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Preface 


This  report  contains  the  results  of  my  efforts  to 
develop  a  computer  code  capable  of  describing  the  time 
dependent  blast  overpressure  and  temperature  profile  result¬ 
ing  from  the  detonation  of  a  nuclear  weapon  in  the  lower 
atmosphere.  The  code  is  not  fully  working  at  this  time, 
due  in  part  to  the  lack  of  valid  opacity  data  at  the  temper¬ 
atures  of  interest.  The  code  is  highly  modular,  and  can  serve 
as  a  starting  point  should  this  work  be  continued  in  the  future. 

I  would  like  to  thank  my  advisor.  Captain  David  Hardin, 
for  his  guidance  and  his  patient,  frequently  repeated,  explan¬ 
ations  throughout  the  research  period.  Most  of  all,  I  would 
like  to  thank  my  wife,  Carol,  for  her  patience,  love  and 
understanding  throughout  the  period  of  this  work. 

Douglas  P .  Wade 
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Abstract 

RADHOT  is  a  one-dimensional  Lagrangian  radiation- 
hydrodynamics  computer  code  developed  from  a  scheme  described 
by  John  Zinn  of  LASL.  The  transport  scheme  described  by 
John  Zinn  and  incorporated  into  RADHOT  requires  just  two 
passes  thru  the  mesh  to  determine  fluxes  and  energy  deposition 
rates.  The  scheme  assumes  local  thermodynamic  equilibrium 
for  each  cell,  and  also  assumes  that  the  thermal  radiation 
flux  contribution  in  a  given  direction  at  any  point  has  one 
of  two  values  depending  upon  whether  a  ray  in  that  direction 
intersects  the  source  region. 

The  RADHOT  code  is  believed  to  be  "debugged,"  but  suffers 
from  a  lack  of  valid  opacity  and  equation  of  state  information. 
Because  of  this,  no  comparisons  can  be  made.  Also,  the  code 
as  currently  written  is  extremely  slow  and  expensive  to  run. 

A  users  guide  and  a  listing  of  the  code  are  provided  as  a 
starting  point  for  future  work. 
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RADHOT:  A  RADIATION 

HYDRODYNAMICS  CODE  FOR 
WEAPON  EFFECTS  CALCULATION 

I  Introduction 

Nuclear  explosives  are  characterized  by  the  very  rapid 
release  of  large  quantities  of  energy  into  a  small  volume. 
Initially,  the  released  energy  is  primarily  in  the  form 
of  thermal  radiation.  As  the  thermal  radiation  streams  away 
from  the  source,  it  is  continously  absorbed  and  re-emitted 
by  the  surrounding  atmosphere.  With  each  absorption/re¬ 
emission,  a  fraction  of  the  thermal  energy  will  be  trans¬ 
formed  into  material  energy.  The  increase  in  material 
energy  is  manifested  as  a  pressure  gradient  between  the 
"absorption"  region  and  the  surrounding,  undisturbed  region. 

This  pressure  differential  results  in  a  shock  wave  being 
initiated  (Ref  1) .  Large  overpressures  associated  with  the 
shock  wave,  together  with  thermal  radiation  fluxes  are  major 
damage  mechanisms  of  a  nuclear  explosion.  An  understanding  of 
the  physics  involved  is  essential  to  producing  a  model  to 
predict  weapon  effects. 

The  difference  equations  of  "shock"  hydrodynamics  are  de¬ 
veloped  by  Richtmeyer  (Ref  2)  and  are  well  understood.  However, 
at  early  times  after  detonation  when  radiation  energy  is  significan 
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and  certainly  prior  to  -10  sec.  (radiative  expansion 
phase) ,  the  hydrodynamic  equations  must  be  modified  to 
include  radiation  transport  (Ref  3)  . 

The  primary  purpose  of  this  thesis  was  to  develop  a 
computer  code  capable  of  predicting  the  weapon  effects  (over¬ 
pressures,  temperatures  and  radiation  fluxes)  resulting  from 
a  1KT  -  1MT  nuclear  explosion.  The  computer  code  was  to 
account  for  radiation  transport  at  early  times,  between 
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~10  sec.  to  10  sec.  after  detonation. 

An  existing  hydrodynamics  code,  HUFF,  was  used  as  a  start¬ 
ing  point  for  this  thesis,  and  RADHOT  incorporates  some  of  HUFF'S 
subroutines  (Ref  4) .  HUFF  utilized  a  1-D  Lagrangian  mesh  to 

solve  shock  hydrodynamic  problems  in  the  absence  of  radiation 
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heating.  At  times  later  than  10  sec.  after  detonation,  HUFF 
is  accurate  to  within  ~5%.  At  earlier  times,  HUFF's  accuracy 
decreases,  and  because  HUFF  lacks  any  provisions  for  radiation 

transport,  it  is  incapable  of  making  any  predictions  prior  to 
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~10  sec.  RADHOT  incorporates  the  I/O  routines  and  the  Lagrangian 
hydro  scheme  from  HUFF,  as  well  as  a  fit  to  the  Nickel-Doan 
Equation  of  State  (EOS)  for  Air  (Ref  5)  . 

The  radiation  transport  scheme  in  RADHOT  was  developed 
by  John  Zinn  of  LASI,  (Ref  6).  In  Zinn's  scheme  radiation 
fluxes  at  the  cell  boundaries  are  determined  at  each  time 
step  by  making  two  passes  through  the  mesh.  The  radiation 
fluxes  are  then  used  to  determine  source  terms  which  are  input 
into  the  hydro  scheme. 
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Detailed  development  of  Zinn's  transport  scheme  together 
with  the  assumptions  and  approximations  behind  it  are  presented 
in  chapter  II  as  is  the  associated  hydro  scheme. 

This  development  comprises  a  major  portion  of  th’s  thesis. 
Chapter  III  discusses  the  timestep  contraints  in  RADFLO 
together  with  the  equation  of  state  and  opacity  approximations. 
The  organization  of  RADFLO  is  discussed  in  chapter  IV,  and 
chapter  V  discusses  the  results  obtained  and  reasons  for  the 
lack  of  correlation  to  data.  Finally,  chapter  VI  contains 
the  conclusions  and  recommendations  of  this  study. 


II  Radiation  Hydrodynamic  P'inite  Difference  Scheme 


Theory 


"Thermal  radiation"  is  electromagnetic  radiation  emitted 
by  matter  in  a  state  of  thermal  excitation.  The  energy  den¬ 
sity  of  such  radiation  is  defined  by  the  Planck  formula: 


B(v,T) 


2hv 3 


(e 


kT 


(1) 


where  v  is  frequency,  T  is  temperature,  B  is  the  Planck  function 
2 

(ergs/cm  ) ,  c  is  the  speed  of  light  and  h  and  k  are  the  Planck 
and  Boltzmann  constants  respectively.  The  importance  of  thermal 
radiation  increases  as  the  temperature  is  raised.  At  moderate 
temperatures,  its  role  is  primarily  one  of  transmitting  energy, 
whereas  at  high  temperatures  the  energy  density  of  the  radiation 
field  itself  becomes  important  as  well.  If  thermal  radiation 
must  be  considered  explicitly  in  a  problem,  the  radiative 
properties  of  matter  must  be  known.  A  determination  of  the 
radiative  properties  usually  requires  an  assumption  that 
"local  thermodynamic  equilibrium"  (LTE)  exists.  Whenever  ther¬ 
mal  radiation  is  considered,  the  medium  that  contains  it  in¬ 
evitably  has  pressure  and  density  gradients  and  the  treatment 
requires  the  use  of  hydrodynamics.  Hydrodynamics  with  explicit 
consideration  of  thermal  radiation  is  called  "radiation  hy¬ 
drodynamics,"  (Ref  7). 

John  Zinn  of  LASL  has  developed  a  computation  a]  method 
for  the  time-dependent  radiation  hydrodynamic  environment 
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encountered  after  the  detonation  of  a  nuclear  weapon  in  t: 
lower  atmosphere.  The  growth  of  a  nuclear  fireball  can  be 
computed  within  certain  limitations  using  the  finite  differ*., 
scheme  that  he  describes.  Zinn's  scheme  is  restricted  to 
problems  with  spherical  symmetry.  The  scheme  requires  just  f. 
passes  through  the  mesh  to  determine  incoming  and  out-going 
fluxes,  and  therefore,  the  scheme  should  be  relatively  fast 
on  a  computer. 

Zinn  makes  a  number  of  simplifying  assumptions  in  the 
scheme,  the  first  of  which  is  the  standard  one  requiring  thnJ 
LTE  exists  in  each  mesh  cell  so  that  local  absorbtion  coeffici 
ents  can  be  expressed  as  functions  of  local  temperature  and 
density.  Zinn  disreguards  radiation  energy  density  and  radic 
pressure  terms  in  the  equation  of  state,  accepting  errors 
at  times  less  than  .10  7  sec.  because  of  the  omission.  Zinn 
ignores  photon  time-of-f light  effects,  justifying  the  omissitr 
because  the  transit  time  of  photons  across  a  cell  is  small  ccn 
pared  to  the  e-folding  time  for  the  exchange  of  energy  between 
the  matter  and  radiation  field.  Finally,  Zinn  starts  his 
scheme  by  placing  all  of  the  energy  of  the  weapon  in  the  first 

_7 

cell  at  t  =  10  sec.  after  detonation  and  restricting  the 
thermal  radiation  from  crossing  the  cell  boundary  until  after 
that  time.  This  prevents  any  preheating  of  the  next  cell. 
Again,  the  principal  effect  of  this  is  to  introduce  a  timing 
\  error  into  the  system. 

Zinn’s  radiation  transport  scheme  assumes  that  the  angu.’ 
distribution  of  the  monochromatic  radiation  intensities  (and 
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ti‘  oforc  the  fluxes)  are  either  1)  approximately  isotropic, 
such  as  in  an  "optically  thick"  region  or  in  an  "optically 
thin"  region  bounded  by  two  "optically  thick"  regions,  or 
2)  can  be  approximated  by  a  function  that  has  one  of  two  values 
such  as  encountered  in  the  cool  transparent  region  outside  of 
an  "optically  thick"  source.  The  transport  scheme  determines 
the  inward-going  and  outward-going  fluxes  at  each  cell  boundary 
for  a  large  number  of  frequency  groups.  By  using  a  large  num¬ 
ber  of  frequency  groups,  a  more  accurate  determination  of  the 
radiative  absorption  coefficient  \i'  can  be  made.  An  accurate 
determination  of  y'  is  essential  to  defining  "optically  thick," 
and  therefore  to  defining  the  source  region.  An  accurate 
determination  of  the  source  region  is  important  to  the  over¬ 
all  accurancy  of  the  scheme. 

The  equation  of  radiative  transfer,  also  called  the 
transport  equation,  is  a  mathematical  statement  of  conserva¬ 
tion.  The  equation  of  radiative  transfer  says 


t 


change  in  transfer  of  attenuation  of 

intensity  =  material  -  radiation  energy 
per  unit  energy  to  (conversion  to 

distance  radiant  material  energy) .  (2) 

energy 


with  the  restrictions  already  noted,  the  equation  of  radiative 
transfer  can  be  stated 
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where 


I  =  radiation  intensity  at  frequency  v 

0  =  angle  between  the  radius  vector  and 

the  direction  of  the  ray  contributing 
to  the  intensity 

3s  =  differential  distance  along  the  ray 
in  the  direction  of  travel 

r  =  distance  from  center  of  symmetry 

\i'  =  absorbtion  coefficient  corrected 
for  stimulated  emmission 

=  Planck  function 
See  Fig  1. 

Now,  Pomraning  (Ref  8:27)  shows  that  the  differential 
equation  of  radiative  transfer  (equation  3)  can  be  replaced 
by  an  equivalent  integral  equation.  The  integral  equation 
or  radiative  transfer  is. 
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(Px)  e 


v'(C)Bv(€)e 
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y'(S)dS 

u'(C')dF.' 
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(4) 


where  dS  and  dE,'  are  elements  of  distance  along  the  ray  con¬ 
necting  P^  and  P2-  See  Fig  2.  In  words,  equation  (4)  says 
that  at  any  point  P2,  the  intensity  contribution  at  that 
point  from  any  other  point  P^  is  equal  to  the  initial  inten¬ 
sity  at  P^  in  the  direction  of  P2,  less  any  attenuation  of 
the  initial  intensity  between  and  P^,  plus  any  re-emmission 


of  radiation  in  the  direction  of  along  the  path  from 

to  less  any  attenuation  of  the  re-emitted  radiation. 

The  rate  at  which  the  material  energy  of  the  system  is 
increased  by  the  absorption  and  re-emission  of  radiation  is, 

3E  \  r 00 

■rr—  )  ,  =  /  U'(  f  I  dio— 4 it  B  )dv  (5 

3t  /rad  J  0  \  J 4tt  v  v/ 


where 
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Ida) 


represents  radiation  energy  lost  to  material 
energy 


and 


4nv  'B 


v 


represents  material  energy  being  given  up  to 
the  radiation  field. 


For  the  transport  scheme  being  developed,  the  flux,  F^,  is 
calculated  instead  of  the  intensity  I  .  Since  the  transport 
scheme  is  restricted  to  problems  with  spherical  symmetry, 
the  two  quantities  are  related  by, 

F  -  f  I  (9)sin0  cos9  dG  (6) 

\J  Jr,  V 


For  problems  with  spherical  symmetry,  equation  (S)  can  be 
expressed  more  simply,  since 
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2 it  sin0  d0 . 
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Let 


COS0  =  0), 


therefore. 


F  =  f  col  da) 

V  /  V 
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Now, 


o)*VI  —  y  I  +  v'B  =  0 
v  v  V 


(9) 


where  wVI^  is  the  net  leakage  out  of  a  volume,  -y'l  repre¬ 
sents  the  loss  of  radiation  due  to  absorption,  and  y'B^  is 
a  source  term.  Then, 


81 
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Integrating  equation  (10)  with  respect  to  solid  angle  yields 
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Therefore,  equation  (5)  can  be  written  as 


3E 

at 
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'rad 


dv 


(13) 


With  the  inclusion  of  equation  (13)  ,  the  equations  of 
radiation  hydrodynamics  are,  the  conservation  of  mass  equation. 


3p  = 

at 


-V- (pv) ; 


(14) 


the  conservation  of  momentum  equation, 

-yP --  =  -V*  ( pvv )  -  VP;  (15) 

and  the  conservation  of  energy  equation, 

3E  °° 

=  -V*(Emv)  -  PV-v  -  J  V-Fv  dv.  (16) 

o 

These  equations,  coupled  with  the  transport  equation,  equation 
(4)  are  the  set  for  which  a  finite  difference  scheme  will 
be  developed  to  solve. 

Equations  13-16  are  a  closed  set,  upon  specification 

of  the  equation  of  state  relations: 

P  =  pressure  =  P(p,E  ),  (17) 

m 

T  =  temperature  =  T(p,E  ), 
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(18) 


the  radiative  absorption  coefficients 


y'  =  y'(p,Em,v)  (19) 

and  the  definition  of  the  flux  and  Planck  functions  (equations 
6  and  1) . 
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Radiation  Transport  Finite  Difference  Formulation 


The  finite  difference  mesh  for  the  category  of  problems 

addressed  in  this  thesis  consists  of  a  series  of  100  concentric 

shells.  Initially  the  thickness  of  each  shell  is  constant, 

though  because  a  Lagrangian  hydro  scheme  is  being  implemented, 

cell  boundaries  will  be  continually  readjusted  during  the 

course  of  the  problem  to  conserve  mass.  Initially,  all  except 

the  first  cell  will  contain  air  at  ambient  temperature  and 

pressure.  The  first  cell  initially  contains  the  mass  of  the 

weapon  (1000  kg)  at  a  temperature  and  pressure  consistent 

with  the  weapon  yield.  See  Fig  3  for  the  mesh  labeling 

-7 

scheme.  Beginning  at  10  sec.  after  detonation,  radiation 

energy  is  allowed  to  propagate  through  the  mesh.  The  equations 

governing  the  radiation  transport  will  be  developed  below. 

At  a  specific  instant  in  time,  the  contribution  to  the 

intensity  I  at  a  point  P„  on  the  surface  of  a  volume  ele- 
v  t  I 

ment  from  any  other  arbitrary  point  on  the  same  surface 
in  a  given  direction  is: 

Iv(P2)  =  Iv(P1)e~y  As  +  Bv(l-e"y  As)  (20) 

as  illustrated  in  Fig  4 . 

This  is  an  approximation  to  the  more  exact  equation  (4) 

r- 

developed  in  the  last  section.  Specifically,  J  p  (£)d£ 

Pi 
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Equation  of  Radiative  Transfer, 


rP2  ^  -f*?  v'  a')dr' 

is  approximated  by  u'is  and  /  B  e  '  d'' 


is  approximated  by  ).  The  dependence  of  u ' 

and  on  the  differential  element  of  length  £  disappears 
because  of  the  assumption  of  local  thermodynamic  equilibrium. 
0^  and  02  are  shown  in  Fig  5,  where  a  ray  is  depicted 


crossing  the  cell.  9^  and  02  are  related  via 


sinO^  =  R2  sin02 


Letting  I1<91)  and  !2*92*  be  the  intensities  of  a  9iven  ray 
at  locations  R^  and  R2  and  As(92)  be  the  distance  between  R1 
and  R2  along  the  ray,  then  an  angle  0m  can  be  defined  for 
b2  such  that 


0  =  sin  (R. / R_ ) 

m  l  z 


which  is  the  maximum  value  of  02  for  rays  that  still  inter- 

o 

sect  R1  (i.e.  the  R1  surface  is  grazed  and  0  ^  =  9  0  )  . 

The  above  equations  can  now  be  used  to  define  the  out¬ 
ward  flux  at  R2 • 

The  out-going  and  in-going  fluxes  at  each  shell  sur¬ 
face  (cell  boundary)  are: 


<  -  2 ”  / 


I  (6)  sin9  cos 9  d9 
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=  —  2 7T  /  I  (G)  sinO  cosO  d9 
J  -n  /  V 
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Using  equations  (18),  (21),  (22),  the  outward  flux  at  R 


-y  'As 


-y  "As. 


Iv(P2)  =  VPl)e  +  V1_e  > 


+  r1 2 

?=2tt  I  (0)  sm9  cos9  d0 

v  J  v 


Therefore, 


=  2t i  |  I^(0^)e  u  As(92)  sj_nQ ^  cos0 


2  d02 


+  27 t  J  B^(l-e  y  As^02^  sin92  cos02  d®2 


+  2n  J  I1(n-02)e  ^  As^02^  sin02  cos02  d0. 
9m 


li/2  > 

+  2tt  J  B^d-e  y  As^e2))  sin02  cos02  d02 


Reducing, 


=  2 7T  J  l^(0^)e  U  As(02)  sj_n0^  cos 


0  dO 

2  2 


B  TTsin  0  -  2tt 


f  m  -y"As  ( 0  2 ) 


sin9_  cosc'  d: 


+  2 


it/  2 

e  ** 

m 

(’f-0.,)  e 

2n 

it/ 

cos  0 

—  2  TT  f 

m 

JQ 

m 

l:'A-S(n.O 


sin-  2  cos' ■  d  ^ 


~U  "  As  ( 0  2 )  . 


sin@2  COSO2  d&2  (26) 


For  the  inward-goi.ng  flux  at  R^  ,  only  those  rays  with 

incident  angles  less  than  0  will  make  a  contribution. 

3  m 

Starting  with 


F.  =  -2tt  f  I  (0  )  cosG.  sinO.  dG  ,  (24) 

1  TT  <  1  1  i. 


and  from  equation  21, 


sin0^  =  r2/,Ri  s^-n®2  (27) 

I(91)  =  I(e2)e~u^As  +  B(l-e~u"As)  (28) 

cos0j  d0^  =  R2^Ri  cos®2  (29) 

2 

sinOj  cos9^  dO^  =  (R2/Rj)  cos02  sinG2  d9^  (30) 


and  from  equation  22, 

(R» /R1 ) 2  =  1 / sin20  (31) 

2  1  m 


therefore , 
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(32) 


sin0^  cos6^ 


sin9»  cos9_/sin  0  dO_ 
2  2  m2 


changing  the  limits  of  integration, 

if  0^  =  tt,  then  ©2  =  tt  (33) 

if  8,  =  tt/2,  then  0  =  tt  —  0  (34) 

J.  2  m 

combining  (24)  with  (28)  ,  and  remembering  that  the  only  rays 
that  make  a  contribution  to  the  flux  at  are  those  whose 
incident  angles  are  less  than  or  equal  to  0^,  then, 

o 

F,  =  -2tt  J  1 0  (  tt  —  6  )e  M  As(02)  sing  cos0_/sin20  d0_ 

1  q  2  2  2  2  m2 

m 

o 

-  2tt  B  J  {1-e  U  As^°2^}sin0„  cos0_/sin20  d0~  (35) 

0  2  4  m  2 

m 

or, 

0m 

F.  =  2/ sin20  f  I  (tt-9  )e  U  As(02)  sinQ  Cos0„  dO 
1  m  J  2  2  222 

o 

9m 

+  ttB  |l-2/sin20m  /  e'lJ  As(°2)  sin02  cos92  d©2  j  (36) 

o 

Now,  it  is  necessary  to  determine  &s(02>.  For  angles  such 

that  0<0_<0  as  shown  in  Fig  6a. 

—  2m 
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03  =  TT-  (TT-0 x )  -  O 2  =  a  -  0 2  .  (37) 

By  the  law  of  cosines, 

As(02)2  =  R2  +  R2  -  2RxR2  cos  (61-02) ,  (38) 

and  using  the  following  identity  (Ref  9) 

cos  ( O 0 2 )  =  cos0^  cos02  +  sinG^  sin@2.  (39) 

Then 

2  2  2 

As(02)  =  R3  +  R2  -  2R3R2  cos0 ^  cos02+sin©1  sin02,  (40) 

2  2 

and  since  sin  0  +  cos  0=1, 

cosO^  =  (1  -  sin201)Jj  (41) 

From  equation  (21)  , 

sin0^  =  R2/R^  (42) 

therefore , 

cos01  =  (1  -  (R2/R:)2  sin202)15,  (43) 
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<>"  f 


cos01  =  1/R^R2  -  R2  sin202)  3s. 
2 

now,  As(02)  is, 

As(02)2  =  _  2R1R2(1/R1 

+  R2/Ri  sin202> 


(44) 


(R2-R2  sin262)J5  cos02 

(45) 


Combining, 

.  2  2 

again,  since  sm  0  +  cos  0  =  1, 

R2  =  R2  (sin202  +  cos202) 

As(02)2  =  r^  +  R2  s^n^®2  +  R2  cos^®2 

-  2R2 (R2  -  R2  sin2@2)  2  cos02  ~  2R2  sin202 
As(02)2  =  R2  -  R2  sin2G2  +  R2  cos 2 ® 2 

-  2R2  cosG  (R2  -  R2  sin202)!5 


(47) 


(48) 


(49) 
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2  2  2  2  2  2  ^ 
As(02)  =  R2  cos  3 2  -  2R2  cos92(R^  -  R2  sin  -j2  )  ^ 

2  2  2 
+  (R^  -  R2  sinz02) 

As(02)  2  =  |r2  cos02  -  (R^-R2  sin202)J5J2 


Therefore 

As ( 0 2 )  =  R2  cos02  -  (R2  -  R2  sin202)^ 

for  O<0_<0 
—  2  m 

For  the  case  where  6m<02_<TT/2 ,  i.e.  a  ray  enters 
at  too  shallow  an  angle  to  enter  R^  before  rc-emerging 
R2,  then  from  Fig  6b,  it  can  be  seen  that 

As(02)  =  21 


and 


i  =  R2  cos02 


Therefore 

As(©2)  =  2R2  COSe2  (0m<02jCTT/2) 


(50) 

(51) 

(52) 

2 

through 

(53) 

(54) 

(55) 
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Zinn  assumes  that  for  each  frequency  group  at  every  point 


outside  the  principal  radiating  sphere,  the  intensity  con¬ 
tribution  for  a  given  ray  direction  has  one  of  two  values. 

If  the  ray  intersects  the  principal  radiating  source  then  the 
intensity  contribution  to  that  point  is  large.  Otherwise  the 
intensity  is  much  less.  From  any  given  point,  the  principal 
radiating  source,  called  the  source  from  here  on,  will  subtend 
angle,  half  of  which  will  be  designated  0  .  Thus  for  points 
exterior  to  the  source,  the  intensity  contributions  at  a  given 
point  will  be  small  for  angles  9  _<6<tt.  See  Fig  7.  The  source 
radius  Rg  must  be  determined  for  each  frequency  group  at  each 
timestep.  Rg  is  defined  as  the  radius  at  which  the  optical 
depth  is  unity,  or 

1  -  zL  ^i(Ri+i  -  V  (56) 

i 


and 


R  =  R 
s  n 


(57) 


where  n  is  the  minimum  number  required  to  meet  the  condition 
in  equation  (56) . 

Once  R  has  been  determined,  0  is  then 
s  s 


9  =  sin"1  (R  /R . ) 

s  si 


(58) 


26 


providing  R^>Rg.  If  R^<Rs,  i.e.  at  a  point  inside  the 

source,  the  intensity  is  then  assumed  to  be  isotropic  and 

0  is  set  equal  to  either  tt/2  or  to  0  =  sin  ^(R.  ,/R  ) 
s  m  l-l  i 

depending  upon  whether  Gg  is  being  determined  as  part  of 
the  calculation  of  an  inward-going  flux. 

The  assumption  that  the  intensity  contribution  at 
each  boundary  has  one  of  two  values  (I  for  0<0  ,  and  I 

3  S  D 

for  9>0g)  permits  the  entire  array  of  fluxes  to  be  generated 
by  just  two  passes  through  the  mesh  (one  inward,  and  one  out¬ 
going  pass)  . 

"This  assumed  intensity  distribution  is  exact  in  two 
limiting  cases;  1)  within  an  optically  thick  region  of 
uniform  temperature  region  of  uniform  temperature, 
where  I  does  not  depend  on  0  (and  then  1^=1^) ,  and 

2)  in  a  vacuum  outside  of  a  uniformly  emitting  surfaces, 
(Ref  6) . 

Now  define 

sin  1  (R  /R. )  for  R. >R 

si  Is 

tt/2  for  R.  <R 

1—  s 

sin  *  (R  /R.  )  for  R.. >R 

si  Is 

0  for  R. <R 

m  1—  s 

Let  and  F+  be  the  inward  and  outward  fluxes  at  R^ 
and  let  F 2  and  F*  be  the  inward  and  outward  fluxes  at  R2 , 

(See  Fig  8)  where  R^  and  R^  are  any  two  adjacent  cell  boundaries 


(59) 


(60) 


x:\  the  entire  mesh  such  that 

Now  using  equations  (23)  and  (24)  and  remembering  that 
I  is  a  constant  with  one  of  two  values  depending  on  9, 


C  * /' 


1,(0)  sinQ  cos9  d9 . 


So,  is 

1  ,v 


F..  =  2tt  I, 

1 ,  v  la 


f  Is  .-TT  /  2 

J  sin9cos9a0  +  2ttI^^  J  sin0  cosG  d0 


or  more  simply, 


Fl,v  =  *Zlasin  9ls  +  7Tllb  (1_sin  0ls)  ' 


and  F_  =  2 tt  I  /  sin0  cos0  d9  +  2ttI_,  /  sin8cos0  d6, 
2'V  29  o  2hJQ2s  (63) 


again  simplifying, 


F2,v  7Tl2a  Sin  J2s  +  12b(1_Sin  62s! 


Starting  from  equation  (24), 


F  =  —  2 tt  /  1.(6)  sinO  cosG  d0 


W 


F,  =  -2i  I...  /  sinO  cos!  dO  =  I., 

i,v  lb  Jv/2  lb 


and 


F_  =  —  2tt  I  .  f  sin6  cos9  d9  =  I„, 

2»v  2b  Jv/2  2b 


[66) 


From  equation  (36) , 


6m 

F,  =  2-r/sin26_  J  I^(7r-0o)e  y  sin6„  cos-3,,  d\ 


m  J  '2  '  2J 
o 


'2  2 


+  ttB(1  -  2/sin29m  j  e  ^  sin02  cosG^  d?2  J  (36) 


F1  “  *lb 


(67) 


Equations  (62)  ,  (64)  ,  (66)  ,  and  (67)  result  from  the  two- 

intensity  assumption,  and  are  simplifications  of  the  more 
exact  equations  (26)  and  (36)  . 

The  integrals  in  equations  (26)  and  (36)  are  all  of  the 

form 


•  a  a  -y'Asi(9)  ,, 

J  s m0  cos 0  e  d0 


(68i 


2  2  .2 

where  As.  (0_)  =  R.,  cosf*.  -  (R.-R0  sin  i^) 

1  Z  Z  Z  X  2* 


(53) 
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or  are  of  the  form 
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33 


-  (R2  -  R1)2/(2R^t2)  ( 1 1  ( :  )  -  H(usp] 

+  7Tllb  |(R1  +  R2)2t2/(4R2)  (G(us)  -  G(um)}J 
-  |(R2  -  R1)2/{2R2t2)  {H(us  -  H(um)}j 

+  Ibb  IR2  -  V2/(2R2l2>  U  -  H'2um'>] 

+  ttB  {1  -  (R  +  R_)2t2/(4R2)  G  (  t  )  -  G(u  )  } 

122  m 

-  |(R2  -  Rj)  2/ (2R2t2)  U  -  H  (  t  )  +  H(2um)}J  (81) 


and. 


F7  =  ttB  +  7T  ( 1 2  b  -  B)  [(Rj  +  R2)2t2/(4R2)  {G(x)-G(um)  } 


-  <R_  -  R  ) 2/ (2R2t2) {H(t)  -  H (u  )} 
2  11  m 


(82) 


Next,  equations  (61)  ,  (62)  ,  and  (63)  can  be  solved  for 
intensities  such  that 


Ilb  F1/tt 


(83) 


I2b  F2/7T 


(84) 
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(85) 


Xla  =  F1  ~  Fl(l  •  sin20ls)/asin2Ols 

I.  =  -  F7/TTsin26.  -  f7/tt  (86) 

la  1  1  Is  1 


Finally,  using  these  equalities,  intensities  can  be 
eliminated  from  equation  (81)  and  (82)  giving  a  result  using 
only  fluxes.  Now, 


F^  =  ttB  +  (F2  -  ttB) 


(R,  +  R0)2t2/(4R2) (G(t)  -  G ( u  ) } 
12  1  m 


(R_  -  R  ) 2/ (2R2t2) (H(t)  -  H (u  ) } 
2  2  2  m 


(87) 


and 

F2  =  +  (Rx  +  R2)2t2/(2R2)  [(F+  -  F”) /sin20ls(G(T)-G(us) } 

+  (F"  -  ttB)  (G(us)  -  G(um)  }J 

+  (R2  -  R1)2/2R2t2)  [(F~  -  ttB)  { 1  -  H(2um)} 

-  (Ft  -  f7) /sin20.  {H (T)  -  H (u  ) } 

1  1  Is  s 

\  -  (F~  -  TIB)  (H(us)  -  H(um)}J  (88) 
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Using  equations  (84)  and  (85)  it  is  now  possible  to 
compute  the  entire  array  of  fluxes  in  just  two  passes. 

The  flux  continutity  conditions  require 


F 


+ 

1 , i+1 ,  v 


2,i,v 


(89) 


and 


F  =  F 

1 , l+l , v  2 , l , v 


(90) 


Starting  at  the  outermost  cell  with  no  inward  flux  (F,  = 

0) ,  equation  (87)  can  be  used  to  compute  all  of  the  inward 
fluxes  in  terms  of  the  previous  flux. 

The  innermost  boundary  of  the  mesh  (r=0)  is  a  point 
in  spherical  geometry,  and  flux  is  thru  an  area,  at  r=0 
both  the  incoming  and  outgoing  flux  is  zero.  With  this  bound¬ 
ary  fixed,  equation  (85)  can  be  used  to  compute  all  of  the 
outward  fluxes  using  quantities  already  computed  during  the 
current  timestep,  for  the  particular  frequency  group. 

Once  the  array  of  fluxes  has  been  generated  for  each 
frequency  group,  the  rate  of  change  of  material  energy  in 
a  particular  cell  due  to  radiation  can  be  determined. 

Recall 


oo 

=  -  J  V • F^dv 

o 


(14) 
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whi  re  _ ml  is  on  a  per  unit  volume  basis.  The  volume  of  any 

Dt  /rad 

particular  cell  is: 


xx 

Ti  =  41T/ 


2  , 

r  dr 


therefore  the  total  rate  of  energy  deposition  in  a  given 
cell  is: 


3Ei  f  o  r  f  *  ? 

_ ml  x  4tt  /  r  dr  =  -  /  V*F  dv  /  4irr  dr  (92) 

3t  /rad  R  \  v  R 


4tt/3  (R^  - 


3E  | 
p3>  — E  ] 

1'  3t/rad 


-/  4*{R2(F2  -  F^l 


-  R^(F+ 


F1) }dv 


Solving  for 


3t  /  rad, 


S^jrad  =  -3/<R2-Rl>/  R2  IP2-P2>  ‘ 


R?<Ft 


F  j )  dv 


Equation  (94)  governs  the  transformation  of  energy  into 
material  energy  which  in  turn  responds  in  a  manner  governed 
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by  the  hydrodynamic  equations.  The  hydrodynamic  scheme  is 

discussed  in  the  next  section,  which  is  followed  by  a  dis¬ 
cussion  of  the  radiation  and  hydrodynamic  timestep  size 
in  chapter  3 . 

Hydrodynamic  Difference  Equations 

The  radiation  transport  method  developed  in  the  last 
section  is  used  to  drive  a  Lagrangian  hydrodynamics  scheme. 

The  radiation  transport  calculations  are  alternated  with  the 
hydro  calculations,  with  one  or  more  hydro  sub-steps  taken 
for  every  radiation  timestep.  The  criteria  for  determining 
the  sizes  of  both  the  radiation  timestep  and  the  hydrodynamic 
timestep  are  presented  in  the  next  chapter. 

Two  hydro  schemes  were  implemented  in  RADHOT.  Originally, 
the  hydro  scheme  was  the  one  used  in  the  HUFF  code.  This 
hydro  scheme  is  a  direct  implementation  of  the  equations  pre¬ 
sented  by  Richtmeyer  and  Morton  (Ref  2:288-324). 

In  this  scheme,  starring  with  known  pressures  and  constant 
initial  mesh  size,  the  size  of  the  hydrodynamic  timestep 
is  computed.  Next,  the  cell  boundary  velocities  and  the  cell 
boundary  positions  are  computed.  Then,  the  cell  densities, 
specific  volumes  and  viscous  pressure  calculated  for  each 
cell.  The  final  step  is  the  calculation  of  the  internal  energy 
and  pressure  of  each  cell.  The  finite  difference  equations 
to  be  solved  in  order  are;  the  cell  boundary  velocity, 


n+1 


u . 


=  u*  +  AtH/p  (P?^ 


n  n 
"  Pi  +  Qj_J 


-  Q *?)/.'  R  <R*?/R0)2,  (95) 
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the  position  of  the  new  cell  boundary. 


_.n+l  _n  ,  n+1  .  , 

Ri  ’  Ri  +  Ui  AtH 


(96) 


the  specific  volume  of  each  cell  (and  hence  the  density) , 


v"+1  -  1/, 


/r>n+1.3  ,  n+1.  3  .  ,„n  .3 

<Ri+1>  -  )  /(Ri+1)  -  (R, 


i»3J 


(97) 


If  the  cell  is  not  expanding  then  the  viscous  pressure  is, 


^n+1  _  2  ,„n+l  ,  _ ,n  ,  n+1  n+1. 

Q.  =  2a  /V.  +  V.  (u.,,  -  u.  ), 

l  i  i  l+l  l 


(98) 


otherwise. 


of1  -  0 


Finally,  in  the  HUFF  hydrodynamic  scheme,  the  next  two 
equations  are  solved  iteratively. 


In+1 

l 


In  - 
l 


P1?*1  +  Pn/2  +  Q*?+1 

l  l  l 


(vn+1 

l 


vn) 

1 


(99) 


P?+1  =  (Y.-l)in+1/Vn+1 
l  '  l  i  l 


(100) 


In  the  equation  above,  is  the  original  cell  density,  Rg 
is  the  original  mesh  spacing,  I  is  internal  energy,  and 
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"a"  is  the  approximate  number  of  cells  over  which  the  shock 
discontinuity  is  spread. 

In  the  process  of  debugging  the  RADHOT  code,  the  above 
hydro  scheme  was  dropped  in  favor  of  one  described  by  Zinn  (Ref 
6)  in  order  to  more  clearly  see  the  coupling  between  the 
radiation  field  and  the  hydrodynamics  it  drives.  The  hydro 
scheme  that  Zinn  describes  for  use  in  conjunction  with  the 
transport  scheme  is  apparently  also  from  Richtmeyer  and 
Morton  (Ref  2)  though  it  is  implemented  in  such  a  manner  as 
to  avoid  the  iterative  procedure  required  by  the  first  scheme 
to  determine  the  internal  energy  and  pressure  of  each  cell. 

The  final  hydrodynamic  finite  difference  equations  and  the 
sequence  that  they  are  to  be  solved  are  discussed  next. 

Initially,  the  pressures,  masses,  densities,  internal 

.  -7 

energy,  and  mesh  spacing  are  determined  for  time  t=10  sec. 

after  detonation.  All  of  the  energy  and  weapon  mass  is  placed 

in  the  first  cell,  with  all  other  cells  undisturbed  from 

ambient  conditions.  The  hydrodynamic  timestep  is  determined 

from  the  Courant  conditions,  and  is  required  to  be  equal 

to  or  less  the  limit  imposed  by  radiation  considerations. 

As  the  time  is  incremented  from  t  to  t+At„,  the  pressure 

ii 

P^  and  the  viscous  pressure  for  each  cell  is  updated. 

The  pressure,  P^,  is  defined  by  the  equation  of  state  and  is, 

p.  =  (y •  -  Dp  i-,  (loi) 

1  '  1  1  1 ' 

and  the  viscious  pressure  is, 
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(102) 


Q.  =  aP.AV./c.At„  {R2  +  R2  } 
1  1  l  1  H  1  1+1 


where  AV^  is  the  change  in  cell  volume,  c  is  the  sound  speed, 
and  a  is  a  constant  between  0.1  and  1.  (fixed  at  .5  in  this 
thesis) .  Equation  (102)  only  applies  to  cells  that  are 
not  expanding.  If  a  cell  is  expanding,  =  0. 


"Next,  the  acceleration  a^+^  of  each  boundary  and  the 
displacement  of  the  boundary  At„  are  computed.  The 
internal  energy  per  unit  mass  for  each  cell,  the  cell 
boundary  location  R^+1  and  the  cell  boundary  velocity, 
,  are  updated.  The  difference  equations  to  be 
applied  in  sequential  order,  are:"  (Ref  6) 


l+l 


=  8-rr  R 


i+1 


(P; 


Qi  - 


i+1 


"W 


(103) 


ARi+l  =  AtH(Vi+l  +  l,ai+lAtH)' 


(104) 


I. 

i 


+  A’<Pi+0i»  lRJ4Ri-RLl  Ri+l>/l"i 


(105) 


i+1 


=  Vi+1  +  3i+lAtH' 


(106) 


and 


Ri+1  =  Ri+1  + 


ARi+l 


(107) 


Then  the  final  cell  volumes  and  densities  are  computed: 
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V. 

1 


4  /  3  7T 


(R 


3 

i+1 


p .  =  m. /V . 

l  li 


A  p  .  =  m .  /  AV . 

i  l  i 


(108) 


(109) 


(110) 


Finally,  the  internal  energies  I.  are  updated  with  At„/At 

1  H 

of  the  radiation  energy  computed  by  equation  (94)  at  the  be¬ 
ginning  of  the  radiation  timestep  At.  Then  prior  to  the  next 
hydro  timestep. 


/3E_ 

I .  =  I .  +  At  !  .  . 
l  l  H\3t 


m 


rad/,4pi»- 


(111) 


at  the  end  of  n  hydro  steps,  all  of  the  energy  in  the 

term  will  have  been  dumped  into  the  hydro  equations 
and  the  next  radiation  timestep  will  be  taken,  beginning  the 
cycle  again. 

After  the  above  scheme  was  implemented,  it  was  checked 
against  the  HUFF  implementation.  The  two  schemes  were  found 
to  agree  to  within  10%.  The  final  scheme  (Zinn's  version) 
should  run  faster  because  it  avoids  the  iteration  required 


in  the  HUFF  scheme. 


Ill  Timestep,  Equation  of  State,  and 


Timestep  Criteria 


In  the  radiation-hydrodynamics  scheme  described  in 
Chapter  II,  there  are  two  different  timesteps  in  concurrent 
use.  The  first  of  these  is  the  radiation  timestep.  At.  The 
rate  at  which  radiation  energy  is  being  lost  to  material  energy 
was  defined  by  equation  (94).  The  radiation  timestep  is  chosen 
such  that  the  fractional  change  in  material  energy  in  any 
cell  in  a  single  timestep  does  not  exceed  101,  or. 


t  =  0.1  min 


SE  \ 

3t  /  i 


(112) 


Zinn  notes  that  this  criteria  is  only  useful  for  problems 
that  are  dominated  by  radiation  that  is  streaming.  If  the 
radiation  transport  is  better  characterized  by  the  diffusion 
approximations,  then  the  timestep  set  by  equation  (112)  becomes 
very  small. 

In  addition  to  the  radiation  timestep,  a  hydrodynamic 
timestep  is  also  used.  The  hydrodynamic  timestep,  At^, 
is  computed  to  conform  to  the  stability  criteria  described 
by  Richtmeyer  and  Morton  (Ref  2) .  Specifically, 


CitH 

1 

Ri+1 


(113) 


where  c^  is  the  local  sound  speed  and  is  defined  (Ref  10) 
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Therefore,  combining  equations  (113)  and  (114),  and  remember¬ 
ing  that  the  criteria  defined  by  equation  (113)  must  be  met 
by  every  cell  in  the  mesh,  then  the  hydrodynamic  timestep 
criteria  becomes, 


R .  .  -  R. 

AttI  <  min  -tV-  1 
H  -  /y.P. 


i  i 
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(115) 
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The  hydrodynamic  timestep  may  not  exceed  the  radiation 
timestep,  although  it  may  be  much  more  limiting.  For  computational 
purposes,  it  is  convenient  to  require  that  the  hydrodynamic 
timestep  be  an  integral  fraction  of  the  radiation  timestep 
so  that  a  fixed  amount  of  the  radiation  energy  can  be  intro¬ 
duced  during  each  of  the  hydrodynamic  substeps.  That  is, 

At  =  At/n,  n=l,2,3,...  (116) 

ri 


To  summarize  the  difference  between  the  two  timesteps, 
the  radiation  timestep  criteria  is  set  primarily  to  insure 
that  the  Planck  function  and  opacity  data  are  updated  frequently 
enough  to  accurately  compute  the  radiation  fluxes  at  each  cell 
boundary,  and  hence  accurately  determine  the  amount  of  material 
energy  being  deposited  by  the  radiation  field.  The  hydrodynamic 
timestep  is  constrained  to  ensure  sufficient  continuity 
of  the  material  properties  of  the  mesh,  and  thereby  avoiding 
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local  instabilities  that  would  destroy  the  finite  difference 
scheme's  convergence. 

Equation  of  State 

The  equation  of  state  for  a  gas  is, 

P  =  (7  -  l)p  I  (117) 

m 

* 

where  P  is  pressure,  p  is  density,  and  I  is  internal  energy. 
Gamma,  y,  is  a  variable  that  ranges  from  1.667  for  an  ideal 
gas  to  a  low  of  approximately  1.05  for  air  at  elevated  temper¬ 
atures.  At  very  high  temperatures  when  a  gas  is  fully  ionized, 
gamma  will  approach  the  ideal  limit.  At  lower  energies,  the 
RADHOT  code  used  a  fit  to  the  Nickel-Doan  Equation  of  State 
for  Air  (Ref  5)  that  was  developed  by  Peters  for  the  HUFF 
code  (Ref  4) .  When  a  significant  fraction  of  the  total  energy 
in  a  cell  is  in  radiation  energy  (>1%)  ,  then  gamma  is  assumed 
to  be  at  the  ideal  limit.  This  approach  turned  out  to  be 
incorrect,  as  it  results  in  a  discontinuity  between  the  regions 
where  Nickel-Doan  applies  and  the  region  where  the  ideal 
gas  approximation  is  valid.  This  discontinuity  can  in  turn 
lead  to  a  splitting  of  the  shock  front. 

Temperature 

The  temperature  of  a  cell  can  be  determined  if  the 
total  energy  in  the  cell  is  known.  Specifically, 

hot  “  I  NkT  +  ir  T* •  <118> 

45 


! '  the  material  energy  I  is  known,  rather  than  the  total 

m 

energy,  then 

T  =  21J  (3Nk)  (119) 

and  a  time  consuming  iterative  procedure  to  determine  cell 
temperature  is  avoided.  In  equation  (119)  ,  k  is  the  Boltzmann 
constant,  and  N  is  the  number  of  particles  in  the  cell. 

For  the  regime  of  temperatures  encompassed  by  the  RADHOT 
code,  N  must  be  able  to  adjust  from  cells  in  which  no  ioni¬ 
zation  has  occurred,  to  cells  in  which  the  air  molecules  are 
fully  ionized,  as  well  as  for  cells  with  all  intermediate 
levels  of  ionization  occurring.  To  do  this  a  fit  was  suggested 
by  Bridgman  (Ref  3)  to  data  reported  by  Zel'dovich  and  Raizer 
(Ref  11) .  This  fit  determines  the  average  number  of  free 
electrons  per  atom,  as  a  function  of  temperature.  See  Fig 
9.  Also  to  be  accounted  for  is  the  disassociation  of  molecular 
nitrogen  and  oxygen  (N^  and  C^)  into  individual  atoms.  There 
are  other  molecular  species  in  the  atmosphere  (NO,  CO,  CO^i 
etc.)  as  well  as  argon,  and  while  some  of  these  species  play 
an  important  role  in  determining  the  opacity  of  the  air,  the 
RADHOT  code  ignors  these  other  species  are  in  the  equation 
of  state.  The  atmosphere  in  the  RADHOT  code  is  assumed  to  be 
79%  N2  and  21%  0^.  The  effective  atomic  weight  of  the  atmosphere, 

A,  is  assumed  to  be  14.59  grams/mole.  Disassociation  of  the 

0 

and  0^  molecules  are  assumed  to  occur  at  20,000  K, 
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v.!'. '  ch  is  also  the  onset  of  where  Z*  (the  number  of  free  electrons 

per  atom)  must  be  included.  Therefore  if  T  is  less  than 
o 

20,000  K,  then, 

N  p 

N  =  — (120) 

In  the  RADHOT  code,  a  variable  Q  is  set  to  either  1  or  2 

depending  upon  whether  the  temperature  is  less  than  or  greater 
o 

than  20,000  K.  From  Fig  9,  Z*  was  determined  to  fit  the 
following  relation: 

4° 

Z*  =  0  ,  T  <  2  x  10  K  (122) 

4  6 0 

Z*  =  4(log10T  -  4.301),  2  x  10  <T<10  K  (123) 

fi  ® 

Z*  =  7  ,  T  >  10  K  (124) 

The  temperature  then  is  defined  as, 

T  =  2ImA/  [3NQPQ(Z*  +  1)]  .  (125) 

Since  Q  and  Z*  are  functions  of  temperature,  correct  deter¬ 
mination  of  each  cell's  temperature  requires  an  iterative 
technique.  This  proved  to  be  too  costly  in  terms  of  computer 
time,  so  Z*  and  Q  are  computed  from  the  old  cell  temperature 
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TEMPERATURE 


first,  and  then  equation  (12:)  is  used  to  compute  the  current 

temperature.  This  method  has  an  inherent  error  associated 
with  it,  but  in  view  of  the  constraints  the  timestep  criteria 
places  on  how  rapidly  the  energy  in  a  cell  is  allowed  to  change, 
the  error  is  considered  acceptable. 

Opacity 

Opacity  is  a  measure  of  the  distance  that  thermal  radia¬ 
tion  will  travel  from  its  source  of  origin  until  it  is  absorbed. 

In  the  lower  atmosphere,  opacity  is  a  complicated  function  with 
minor  atmospheric  constituents  playing  important  roles.  Very 
large  computer  codes  (DANA,  DIANE,  DIAPHANOUS)  have  been  developed 
for  the  accurate  prediction  of  absorbtion  coefficients.  These 
codes  are  maintained  by  LASL  and  the  Air  Force  Weapons  Labora¬ 
tory  (AFWL).  Opacity  data  and  equation  of  state  data  were  r 
quested,  though  never  obtained,  from  AFWL. 

Failing  to  obtain  data  from  AFWL,  a  fit  was  made  to  a 
graph  of  absorption  coefficients  that  appeared  in  Zinn's 
paper  (Ref  6) .  See  Fig  10.  The  values  plotted  in  Fig  10 
are  for  ambient  air  density  only.  When  the  air  density  varies 
from  ambient,  values  obtained  from  fits  to  Fig  10  are  subject 
to  large  (orders  of  magnitude)  errors.  Figures  11  and  12 
(Ref  12)  illustrate  this  point. 

Using  fits  to  Fig  10  data  proved  unsatisfactory  in 
RADHOT  (negative  fluxes  occasionally  resulted) ,  so  the  absorb¬ 
tion  coefficients  were  averaged  over  frequency,  and  a  smoothly 


49 


Fig  10.  Plots  of  radiative 

absorbticn  coefficients 
vs.  photon  energy  (Ref  6) 


Fig.  12 Absorption  (’n«  flu  icnt  <  f  Air  :is  a  Function  of  the  Photon  Knorgj 
v  IP'-*  gir/cnv*;  Trmivr.'Kurr.  -JlhifK  (Ref  12) 


varying  function  of  temperature  was  obtained.  Specifically 
the  Rosseland  mean  is  being  computed,  and  it  is. 


(126) 


where  B  is  the  Planck  function  (equation  (1)).  After 

3° 

y  was  determined  for  various  temperatures  between  10  K 
7°  - 

and  10  K,  then  y  can  be  determined  for  any  temperature  in 


that  range  by  the  following  relation: 


An[jj'(T)  1]  =  4n[p' 
+  *n  ^  |(Ti+l)_1) 


(Ti) 


-1 


in  Ti+1  -  in  T 

in  T.  ,  ,  -  inT." 
l+l  l 


inT  -  in  T. 

_ l 
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(127) 


The  Rosseland-averaged  absorption  coefficient  when  used 
in  RADHOT  no  longer  predicts  negative  fluxes,  but  since  the 
\i'  used  in  determining  the  averages  was  from  Fig  10,  variations 
in  density  can  still  result  in  calculated  values  of  p'  which 
may  be  off  by  orders  of  magnitude.  This  second  approach  results 
in  the  absorption  coeffieient  being  strictly  a  function  of  temp¬ 
erature.  At  the  high  temperatures  that  are  encountered  in  the 
innermost  cells  at  the  start  of  the  problem,  the  absorption 
coefficient  forces  energy  transport  to  be  of  a  diffusive  nature 
(the  cells  are  optically  "thick").  The  Rosseland  averaging 
effectively  results  in  the  existance  of  a  single  frequency 
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nr  an.  With  a  single  frequency  group  there  is  no  possibility, 
through  absorption  and  re-radiation,  of  a  portion  of  the  energy 
in  the  cell  being  re-emitted  at  a  frequency  where  the  cell  is 
not  optically  "thick"  and  hence  continuing  the  streaming  pro¬ 
blem.  Further  discussion  of  the  impact  of  the  lack  of  valid 
opacity  data  on  the  RADHOT  code  is  in  chapter  V. 
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IV  Problem  bef in i ti on  and  Sol utlon 

Program  RADHOT  is  designed  to  solve  one-dimensional 
spherical  geometry  problems  associated  with  the  detonation 
of  a  nuclear  weapon  in  the  lower  atmosphere.  The  program 
describes  the  deposition  of  thermal  radiation  energy  from  the 
weapon  into  the  surrounding  air  resulting  in  the  development 
and  propagation  of  a  shock  front.  The  program  tracks  the  temper¬ 
ature,  pressure,  density,  boundary  velocity,  and  viscosity  of 
each  of  the  100  mesh  cells  as  a  function  of  time,  printing  out 
the  information  for  times  requested  by  the  user. 

Program  RADHOT  retains  the  modular  programming  concepts 
of  its  predecessor,  HUFF,  as  well  as  several  subroutines  from 
the  original  code.  The  retained  subroutines  have,  in  most 
instances  been  modified  for  use  in  RADHOT.  RADHOT  is  written 
in  FORTRAN  IV  as  implemented  on  the  CDC  CYBER  series  computer. 

The  radiation  hydrodynamics  problem  is  defined  by  sub¬ 
routines  BLAST  and  RADSET.  The  initial  conditions  are  based 
on  the  following: 

1.  The  weapon  yield  is  deposited  in  the  first  mesh  cell, 
and  is  divided  between  a  thermal  radiation  field,  and  the  in¬ 
ternal  energy  deposited  in  the  mass  of  the  first  cell. 

2.  The  mass  in  the  first  cell  is  the  mass  of  the  weapon, 
and  is  arbitrarily  fixed  at  1000  kg  (10^  grams) .  The  remain¬ 
ing  99  cells  in  the  mesh  contain  ambient  air  at  a  pressure 

and  density  adjusted  for  altitude. 

3.  The  first  cell  is  isothermal. 


55 


-7 

4.  The  problem  begins  at  10  seconds  after  detonation 
(assumed  to  be  the  end  of  the  neutronics  phase) . 

Subroutine  BLAST  is  responsible  for  reading  in  the  weapon 
yield  in  KT  and  altitude  scaling  factors  for  different  height 
or  burst  (HOB) .  The  subroutine  scales  the  programmed  values 
calculating  the  proper  initial  conditions  for  a  variety  of 
problems.  The  altitude  scaling  parameters  provided  in  Table  I 
are  input  by  the  user.  The  parameters  are  defined  as 


(128) 


(129) 


(130) 


2 

where  PQ  and  TQ  are  1.013  dyne-cm  and  300 °K  respectively  (Ref 
1:104).  Subroutine  BLAST  creates  the  100  cell  mesh,  defining 
initial  cell  boundaries,  densities,  pressures,  and  cell  volumes. 

Subroutine  RADSET  is  responsible  for  further  initializ¬ 
ing  values  in  the  mesh  for  the  radiation  transport  problem. 
Subroutine  RADSET  reads  the  number  of  frequency  groups  to  be 
used  and,  for  every  cell  except  the  first  one,  RADSET  determines 
the  mass  of  air  in  the  cell,  and  fixes  the  amount  cf  ambient 
radiation  and  material  energy  in  the  cell  on  a  per  volume 
basis  to  a  level  consistent  with  a  temperature  of  300°K. 

-7 

RADSET  also  fixes  the  start  time  of  the  problem  to  10 

sec.  For  the  first  mesh  cell,  RADSET  in  turn  calls  another 
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Table  I 

Average  Atmospheric  Paramters  (Ref  1:104) 


Altitude 

(feet) 

Scaling  Factors 

Altitude 

(feet) 

Scaling  Factors 

SP 

SD 

ST 

SP 

SD 

ST 

00 

1.00 

1.00 

1.00 

• 

15,000 

0.56 

1.21 

1.28 

1,000 

0.9C 

1.01 

1.02 

20,000 

0.46 

1.30 

1.39 

2,000 

0.93 

1.03 

1.03 

25,00' 

0.37 

1.39 

1.53 

3,000 

0.90 

1.04 

1.05 

30,000 

0.30 

1.50 

1.68 

4,000 

0.05 

1.05 

1.07 

35,000 

0.24 

1.62 

1.86 

5,000 

C.83 

1.06 

1.08 

40,000 

0.19 

1.75 

2.02 

10,000 

0.69 

1.13 

1.17 

57 


'outine,  subroutine  TEMPO,  to  determine  the  temperature  and 

the  amount  of  radiation  and  material  energy  in  the  first  cell. 
Currently,  RADFLO's  initializing  subroutines  are  fixed  to 
provide  information  for  a  one  megaton  (1MT)  weapon,  with  a 
weapon  mass  fixed  at  1000  kg.  A  1MT  weapon  results  in  an 
initial  mesh  spacing  of  89  cm.  Modification  to  the  initial¬ 
ing  subroutines  (BLAST  and  RADSET)  would  be  required  for  pro¬ 
blems  with  different  yields  and  weapon  masses. 

The  RADHOT  program  is  depicted  in  Fig  13.  The  program 
begins  by  calling  subroutines  FLAGS.  Subroutine  FLAGS  deter¬ 
mines  at  what  times  output  is  requested.  RADHOT  then  trans¬ 
fers  control  to  subroutine  CONTROL. 

Subroutine  CONTROL  is  responsible  for  initializing  the 
program,  directing  the  flow  of  the  program  and  printing 
the  results.  Subroutine  CONTROL  first  sets  the  initial  con¬ 
ditions  by  calling  BLAST  and  RADSET  and  sets  up  a  loop  in  which 
the  radiation  hydrodynamic  calculations  are  accomplished. 

First,  subroutine  CONTROL  calls  subroutine  Planck  which 
evaluates  the  Planck  function  for  each  cell  at  each  frequency. 
Then,  using  opacity  data  generated  by  function  AMU,  the  size 
of  the  principal  radiating  source  at  each  frequency  is 
determined.  These  last  two  steps  generate  input  that  will 
be  required  for  the  radiation  transport  calculations.  Next, 
the  radiation  and  hydrodynamic  timesteps  are  determined, 
using  criteria  presented  in  chapter  III. 

Then,  subroutine  RADFLO  is  called.  Subroutine  RADFLO 
is  direct  implementation  of  the  radiation  transport  equations 
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If  cycle  I 

limit  not  Determine  At 

reached  I 


Fig  13. 


Program  RADHOT  Flow  Diagram 


developed  in  chapter  II.  Subroutine  RADFLO  determines  the 
net  amount  of  radiation  energy  being  dumped  into  material 
energy  in  each  cell  during  the  current  radiation  timestep. 
After  RADFLO  is  executed,  subroutine  HYDRO  is  executed  one  or 
more  times.  Subroutine  HYDRO  implements  the  second  hydro 
scheme  described  in  chapter  II.  At  the  end  of  each  hydro- 
dynamic  substep,  the  internal  energy  in  each  cell  is  updated 
with  the  appropriate  fraction  of  the  radiation  energy  being 
deposited,  and  the  pressure  and  temperature  of  each  cell 
is  recomputed.  Cell  temperature  is  computed  in  a  separate 
subroutine  named  TEMPER.  Subroutine  TEMPER  incorporates  the 
Z*  approximation  described  in  chapter  III  via  function 
ZSTAR . 

At  the  end  of  each  radiation  timestep,  the  elapsed 
time  is  updated,  and  the  location  of  the  shock  front  is 
identified.  If  the  shock  front  has  reached  the  94th  cell, 
then  subroutine  REZONE  is  executed  and  the  cell  boundaries 
are  relocated  in  such  a  manner  as  to  conserve  energy. 

The  cell  boundaries  are  moved  such  that  the  shock  front 
after  re-location  is  in  cell  49. 

Finally,  prior  to  continuing  on  to  the  next  radiation 
step,  a  check  is  made  to  determine  if  it  is  time  for  printing 
results.  If  it  is,  then  subroutine  WRITER  prints  cell  para¬ 
meters  (boundaries,  pressures,  temperatures,  densities,  radia¬ 
tion  and  internal  energy  in  the  cell,  and  in-going  and  out¬ 
going  fluxes  at  the  cell  boundaries.  Line  printer  riots 
of  cell  parameters  vs  position  is  accomplished  fcv  subroutines 
COPY  and  PLOTSM. 
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The  program  continues  cycling  until  either  the  time  limit 
or  the  cycle  limit  specified  by  the  user  (in  subroutine  FLAGS) 

is  reached. 
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V  Res u Its  and  D.i  scussi  on 

Program  RADHOT  is  believed  to  be  "debugged"  but  is  • 
producing  reasonable  results  due  to  a  lack  of  good,  freer, 
dependent,  opacity  data,  and  due  to  an  incomplete  equatir: 
state.  The  equation  of  state  and  opacity  data  currently 
in  the  code  shows  a  shock  front  developing  and  propagating 
orders  of  magnitude  slower  than  expected.  The  subroutine  J 
calculates  the  absorption  coefficient  is  predicting  values 
in  the  vicinity  of  the  innermost  cells  that  are  far  too  s • 
The  result  is  the  RADHOT  code,  which  is  optimized  to  hanc . 
streaming,  is  calculating  fluxes  in  a  diffusive  region. 

RADHOT  has  provisions  for  treating  optically  thick 
cells  with  diffusion  theory,  but  the  code  expects  that  me. 
of  the  energy  transfer  thru  a  cell  will  be  occuring  at 
frequencies  where  the  cell  is  not  optically  thick.  Unfor' 
unately,  the  averaging  of  the  absorption  coefficients  des¬ 
cribed  in  chapter  III  removed  the  frequency  dependence  and 
made  the  absorption  coefficient  strictly  a  function  of  tem¬ 
perature.  The  result  is  that  once  an  optically  thick  cell 
is  reached,  there  is  no  way  through  absorbtion  and  re-emis. 
at  a  different  frequency  to  stream  through  that  cel]  .  The 
energy  transfer  occurring  is  by  diffusion.  Zinn,  in  a  com- 
similar  to  RADHOT,  used  40  frequency  groups.  C).  R ' 

has  provisions  in  it  for  up  to  eight  frequency  gr:i.:.s,  1,  i.  _ 
to  the  current  state  of  the  opacity  subroutine,  there  is 
advantage  to  running  more  than  one  group. 

Another  problem  that  arises  when  the  code  is  in  a  di 
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regime  is  that  the  radiation  timestep  criteria,  equation 
(112),  computes  timesteps  that  are  inconviently  small.- 
During  the  development  of  the  program,  the  timestep  was  noted 

to  grow  increasingly  smaller  as  the  program  progressed. 

-20 

If  allowed  to  run  long  enough,  timesteps  of  10  sec.  were 
eventually  generated. 

In  order  to  verify  that  the  transport  scheme  developed 
in  chapter  II  was  working,  the  portion  of  the  code  responsible 
for  diffusion  theory  calculations  was  "switched  off,"  and 
the  timestep  size  was  fixed.  Doing  this  shows  that  the  trans¬ 
port  code  can  calculate  in  a  diffusive  regime.  The  shock 
front  slowly  propagates  outward,  and  the  mesh  values  remains 
stable  and  physically  "reasonable,"  (no  negative  temperatures, 
energies,  or  fluxes  and  all  quantities  move  in  the  expected 
directions  as  time  progresses) .  Figure  14  is  the  printout 
of  computed  parameters  of  the  first  20  cells  for  the  first  two 
cycles.  By  the  second  cycle  the  second  mesh  cell  is  beginning 
to  show  the  increase  in  material  energy  as  manifested  by  a 
slight  increase  in  cell  temperature  and  pressure.  See  Fig 
14. 

At  7.6  x  10  4  sec.  after  detonation  (850  cycles)  the  rise 
in  temperature  and  pressure  in  the  second  cell  are  more  pro¬ 
nounced,  and  the  DEDT  column  indicates  that  energy  is  beginning 
to  be  dumped  into  the  third  cell,  though  the  effect  is  not  yet 
noticeable.  See  Fig  15.  Also  at  this  point  in  time,  the  DEDT 
term  shows  a  very  small  amount  of  energy  is  bei"  i  deposited 
in  cells  5  and  6.  This  pair  of  DEDT  terms  detach;-.-  Itself 
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f>'  the  principle  radiating  source  and  travels  rapidly  through 
the  mesh.  This  phenomena  occurs  at  random  intervals  during 
the  program  run.  This  may  be  the  "shock  splitting"  mentioned 
in  chapter  III,  arising  due  to  the  discontinuity  in  the  equa¬ 
tion  of  state,  or  it  may  merely  be  a  computational  anomaly. 

By  2.7  x  10  ^  sec.  after  detonation  (cycle  2800),  the 
first  two  cells  are  beginning  to  come  into  equilibrium  with 
each  other.  Also  the  boundary  between  the  first  and  second 
cell  is  beginning  to  expand  outward,  compressing  the  air  in 
the  second  cell.  See  Fig  16.  Finally,  by  4.5  x  10  sec. 
(cycle  460 0)  the  first  two  cells  are  nearly  in  equilibrium  with 
each  other  and  the  third  cell  continues  to  heat  up.  See  Fig 
17.  The  reason,  in  part,  for  the  slowing  of  the  shock  front 
is  that  though  each  concentric  shell  has  approximately  the  same 
thickness,  the  total  amount  of  mass  of  air  in  each  shell  in¬ 
creases  with  radial  distance  from  the  source. 

Due  to  the  amount  of  CPU  time  required  to  run  RADHOT, 
the  problem  has  not  been  run  beyond  4600  cycles.  Also,  the 
above  run  was  for  a  single  frequency  group.  RADHOT  required 
800  seconds  of  CP  time  to  run  4600  cycles,  or  roughly  0.18 
seconds  per  cycle.  Running  with  more  than  one  frequency  group, 
as  would  be  required  to  do  the  problem  properly,  would 
further  worsen  the  situation.  Further  work  on  RADHOT  will, 
in  addition  to  improving  the  opacity  data,  have  to  look  at 
ways  of  decreasing  the  cycle  time.  RADHOT  does  not,  however 
require  large  amounts  of  memory,  approximately  68,000  memory 
locations. 
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Finally,  one  further  check  was  made  to  check  the  function- 

-3  -1 

mg  of  the  program.  The  opacity  was  fixed  at  10  cm  through¬ 
out  the  run.  The  results  are  shown  in  figures  18  and  19. 

By  the  second  cycle,  energy  is  streaming  strongly  outward  as 

denoted  by  the  F(+)  column.  The  outward  flux  falls  off  roughly 
2 

as  1/R  .  By  cycle  850  energy  is  still  flowing  through  the 
mesh  as  denoted  by  the  DEDT  terms.  The  temperature  of  the 
innermost  cell  is  lower  than  in  the  previous  example,  and 
this  is  due  to  the  energy  that  was  lost  when  it  streamed  to 
be  deposited  in  cells  further  out.  This  result  is  what  would 
be  expected  if  the  code  was  running  properly,  but  contained 
bad  opacity  and  equation  of  state  information. 


Fig  19.  1MT  BLAST,  FIXED  p',  1.  6  x  10  sec 


Conclusion 


VI  Conclusions  and  Recommendations 


The  RADHOT  program  implements  a  radiation  transport 
scheme  coupled  with  a  hydrodynamics  scheme  that  is  reported 
to  be  both  fast  and  fairly  accurate  in  calculating  the  en¬ 
vironment  around  a  nuclear  weapon  at  early  times.  The 
RADHOT  code  is  currently  unable  to  evaluate  the  preceding 
statement  due  to  a  lack  of  good  opacity  data.  The  opacity 
data  currently  in  RADHOT  makes  the  air  around  the  source 
region  highly  opaque,  resulting  in  a  code  optimized  for  sit¬ 
uations  characterized  primarily  for  streaming,  attempting 
to  calculate  the  behavior  of  a  diffusive  regime.  Correction 
of  the  opacity  data  is  essential  to  any  evaluation  of  the 
accuracy  of  the  RADHOT  code.  The  equation  of  state  used  in 
RADHOT  also  needs  further  work,  i.e.  the  "Nickel-Doan  Equation 

r  0 

of  State  for  Air"  fit  needs  to  be  extended  to  10°  K. 

The  original  approach  in  this  thesis  was  to  attempt  to 
couple  the  radiation  transport  scheme  described  by  Zinn  to 
an  existing  hydrodynamics  code  (HUFF) .  A  great  amount  of 
effort  and  time  was  required  to  understand  and  debug  the  HUFF 
code.  Further,  the  final  code  resulting  from  this  "forced 
marriage"  is  inefficient,  for  example,  in  the  RADHOT  code  there 
are  arrays  of  cell  volumes  and  mass,  cell  densities,  and 
cell  specific  volumes  because  the  different  parts  of  the  program 
require  different  forms  of  the  same  quantity.  Computing 
redundant  quantities  partially  accounts  for  the  slow  cycle 
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t  -  in  RADHOT.  Using  a  look-up  table  for  the  Planck  function, 
as  well  as  for  the  opacity  data  would  also  accelerate  the 
cycle  time. 

A  better  approach  would  have  been  to  program  Zinn's 
radiation  hydrodynamics  scheme  directly,  allow  it  to  run  until 
shortly  after  hydroseparation  (the  starting  point  of  HUFF) , 
and  then  provide  output  compatible  with  the  input  requirements 
of  HUFF.  By  running  the  two  programs  in  series,  the  radiation 
hydrodynamics  could  be  turned  off  when  it  is  no  longer  signi¬ 
ficant,  and  HUFF  would  have  its  initial  conditions  corrected 
for  radiation. 

Recommendations 

1.  Obtain  valid  opacity  data  for  RADHOT.  This  is  the 
single  most  important  item  that  needs  to  be  accomplished. 

2.  Improve  the  equation  of  state  fit  by  smoothly  extend¬ 
ing  the  Nickel-Doan  fit  of  gamma  to  the  ideal  limit  of  5/3. 

3.  Modify  RADHOT  so  that  at  the  completion  of  a  run, 
it  will  provide  input  data  to  HUFF. 

4.  Incorporate  look-up  tables  for  the  Planck  function 
and  for  opacity  data  into  RADHOT. 

5.  Incorporate  programing  techniques  into  RADHOT  that 
can  increase  the  execution  speed. 
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Appendix  A 


Derivations 


As(0)  =  R2  cosG 


AS(02s)  =  ? 


?  2  2  k 

(R^  -  R2  sin  a) ^ 


92s  =  Sin_1  f°r  R1  '  Rs 


2  2  k 

As(02s)  =  r2  cos02s  ~  (R1  “  V 


As  ( 0  2g) 2  =  r2(1  -  sin^02s)  "  R2  COS°2s  (R1  "Rs)  '  +  Rl~Rs  (A-4) 


(A-l) 


(A-2) 


(A—  3 ) 


.2  _2 .  *5  .  „2  „2 


As  (0  ) 

2s 


n  rs. 

R2  (1-  -§) 
R2 


2  y  2  2 

R  cosO,  (R,-R  )  2  +  (R-.-R  ) 
2  2s  1  s  is 


( A—  5 ) 


4s<e2s)2  =  (R2-R2)  -  R2  COS02S(R2-R2)’i  + 


isle  j2  =  (Rj-rJ)  -  iR2  cos2°2stRrI,2,,,i  +  (RrRi* 


4S(02S) 


*s(02s> 


2  2 
<R2-Rs> 


2  2 


(R^(l-sin202s) (R2-R2) )h  +  (Ri-Rg) 


(r2-rJ)  h  +  (Ri'R^ 


'1  “s' 


:4S(0  )  »  <<r2-r2>  -  (R?-Rj))!i 


u(02s)  =  u'as(02s>  =  -fR^qy  4sI02Si 


(A-6) 


(A-7 ) 


(A-8 ) 


(A-9) 


(A-10) 


(A- 1 1 ) 
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(73) 


u  (  'o  J  '  ( 


:/AR)((R^-R^  "  (Ri  R2)  J 


u  =  u(6_) 
m  m 


.  -1  R1 
q  =  sm  r~ 
m  k2 


(A-12) 


u  ( 0 )  =  U'AS(0)  =  v.MR2  cose  -  <r\  4  sin2e)2) 


(A-13) 


As(Qm) 


R2  cos9m  * 


2  2  .  2 .  >  h 

(Rj~R2  sin  V 


(A-14) 


As(em) 


As  (0  ) 
m 


R2  cos8m  - 


R2  COS0m 


>r?-R2  ^  )N 

«RrRi> 


(A-15) 


(A-16) 


As  ( 0  )  =  R2  cosGt 


(A-17) 


As(0m) 


cos2 e  -  R2  l1 
2  m  z 


sin2°m’ 


(A-18) 


As|V 


r4  » 


(A-19) 


9  2  2 
As(0)  =  r2  ~  R1 


(A-20) 


As  (0, 


)  =  Ur2  +  Rl)  (R2  '  Rl}) 


(A-21) 


u,0  )  =  U'ASHJ  *  (x/(82-Ri"<<R2+RlMVRl,) 


(A- 22) 


u ( 0  )  =  T < <R2+R1) / (R2'R1) * 
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Appendix  B 

User's  Guide  for  RADHOT 


This  appendix  describes  the  data  card  input  and  the 
NOS/BE  control  cards,  and  selecting  the  output  when  operating 
program  RADHOT.  The  data  card  input  will  be  presented  first 
followed  by  a  discussion  of  the  NOS/BE  control  cards  required 
to  run  RADHOT. 

Data  Cards  for  RADHOT 

RADHOT  is  designed  to  solve  the  problem  of  a  nuclear 
detonation  in  air.  Due  to  the  initial  mesh  spacing  being 
fixed  in  the  program  (subroutine  BLAST) ,  the  program  as  it 
currently  exists  can  only  solve  this  problem  for  1MT  of 
yield  at  sea  level.  Addition  of  a  scaling  parameter  for  the 
initial  mesh  spacing  will  allow  the  program  to  be  run  for  dif¬ 
ferent  yields  and  altitudes.  This  users  guide  will  define  what 
data  is  required,  and  what  the  value  of  the  data  is  currently 
fixed  at. 

Parameters 

1.  Card  1-FLAG  (2),  FLAG (3),  FLAG (4),  TSTOP 

This  card  contains  the  input  flags  defined  as  follows: 
FLAG (2)  is  an  integer  representing  the  number  of 
result  times  listed  on  card  2. 

FLAG  (3)  is  an  integer  representing  the  rraxium  number 
of  timestep  iterations.  The  program  will  stop  if 
it  reaches  this  many  iterations. 

FLAGS (4)  is  an  integer  representing  the  type  of  problem 
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to  be  solved.  FLAG (4)  must  be  3. 

TSTOP  is  the  maxium  problem  time  allowed.  The  program 
will  stop  if  it  reaches  this  time.  Results  for  this 
time  will  be  plotted  in  addition  to  those  on  card  2. 

2.  Card  2-TABLE (1),  TABLE ( 2 )  ,  .  .  .  ,  TABLE (FLAG (2) )  .  This 
card  lists  the  requested  result  times.  A  maxium  of  10 
times  can  be  listed  due  to  the  dimension  of  the  array 
TABLE . 

3.  Card  4-MAXNU.  Number  of  frequency  groups  used  in  the 
problem.  Maxium  value  is  8,  though  because  of  the  current 
poor  opacity  data  in  subroutine  AMU,  there  is  no 
advantage  to  having  more  than  one  group. 

Problem  Defining  Card 

Card  3-YKT ,  SP,  DS,  ST  is  the  problem  defining  card. 

The  variables  are:  YKT-Yield  of  the  blast  in  kilotons. 

currently  fixed  at  1000  ( 1MT) 

SP,  DS,  ST-altitude  scaling  factors  from  Table  1  of  the 

text.  Currently  fixed  for  sea  level. 

c 

Current  Data  Set 

The  following  data  cards  are  being  run  with  RADHOT: 

5  10000  3  1. 

l.E-6  5.E-5  3.E-4  9.E-2  .8 

1000.  1.  1.  1. 

1 

NOS/BE  Control  Cards  for  RADHOT 

This  section  presents  the  control  cards  necessary  to 
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execute  the  RADHOT  program  at  AFIT. 


Permanent  File.  To  place  RADHOT  on  permanent  file,  the 

following  cards  are  necessary: 

job.  problem 
REQUEST,  RADHOT,  *PF. 

COPYBR ,  INPUT,  RADHOT. 

CATALOG,  RADHOT,  RP=999. 

7/8/9 

RADHOT  deck. 

7/8/9 
6/ 7/8/ 9 

after  RADHOT  is  on  permanent  file,  the  following  control  cards 

are  used  to  run  the  program. 

job,  T200,  CM75000.  problem 
ATTACH, IMSL,  ID=LIBRARY ,  SN=ASD. 

LIBRARY, IMSL 
ATTACH, RADHOT. 

FTN , I=RADHOT , 0PT=2 . 

LGO . 

7/8/9 

DATA  cards  for  RADHOT 

7/8/9 

6/7/8/9 . 
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Appendix  C 


1. 00 
Ho 
1 2  0 = 
130: 
140= 
150=: 
160--: 
170:= 
180  = 
190  C 
200 
210= 
220= 
230=0- 
240  = 
250  = 
260= 
270  = 
200= 
290=-  C 
300=C 
310==C 
320--C 
330== 
340“ C- 
350== 
360= 
370==C- 
380== 
390  = 
400  = 
410  = 
420== 
430  = 
440  = 
450= 
460= 
470== 
480  = 
490== 
500  = 
510  = 

[  520== 

I  530  = 
540  = 


RAD HOT  Program  Listing 


.i  i I...L5 > review ,■  ;  j lls >■  fori:;,  review; 

H  I  i  1 1  N  8  i:  0  N  I  L  A  0  <  5 )  »  T  A  B  I,  IT.  (  J.  0  ) 

INTEGER  FLAG 

CALL  FLAGS ( I  I  AG  y I ABLE  y TSTOP ) 

C  A  I...  I.  C  G  N  T  R  D  I...  ( i:  I  A  G  y  1  A  l’-l  l.y  T  B  T  0  P  ) 

PR  I  NT  if  "I  I  ME  L.  I M I T  EXCEEDED  “ 

STOP 

END 


S U B ROUT  I N E  F I .  A G S  <  F L  .A G  v  T  ABLE »  TS  T  0  P  ) 

DIMENSION  FLAG ( 5 ) y TABLE < 10 ) 

INTEGER  FLAG 

-  INITIALIZE  CONTROL  PARAMETERS  - 

DO  10  1=1 y 3 
10  FLAG  <  I. )  =  0 
FLAG (5)  =  1 
DO  30  1=1 » 10 
30  TABLE ( I ) =0 

I  READ  FLAGS  AND  REAL  TIME  LIMIT  I 

I  FLAG  ( 2 )  =  11  OF  TABLES  REQUESTED  I 

I  FLAG ( 3 ) ==#  OF  CYCLES  THROUGH  HYDRO  I 

I  FLAG  (  5 )  =  1  v  INDICATES  l-IRST  PROGRAM  ITERATION 

READY y ( FLAG <  I )  » 1  =  2 y  4 ) y TSTOP 

-  READ  REOUESTED  TABLE  TIMES  - 

IF 2- FLAG (2) 

READY y  <  TABLE ( I ) y I  =  1 v I F2 ) 

- WRITE  OUT  INFORMATION  PROVIDED  BY  FLAGS  - 

WRITE ( 2  y  400 ) 

WRITE ( 2  y  700 ) 

W  R I TE  <  2  y  200 )  F I ...  A  G  (  2 )  v  ( T  A  B I...  E<I)»I  =  1»IF2) 

W R I T E <  2 y 3 0 0 ) T S f 0 P  y  FLA G ( 3 ) 

2 0 0  F 0 R M A T  (////»!’  1 0  y  “  T A B L E S  F 0 R  T I  I E  I'" 0 L. L. 0 U I N G  ’  ?  1 3  r 

1 ‘  T I M E S  W E R E  R E Q U E S T E D I “ y // y  2 ( T 1 5 y  5 ( 1 P E 1 0 . 3 » 2 X ) y / ) ) 

300  F  0  R  M  A  T  <////»  T 1 0  y  "  M  A  X  T I M  E  T  0  R  LJ  N  I.  S  “  ,  F  4 . 1  y  “  SEC." 

1 y/yT10y "Max  CYCLES  THROUGH  TI HESTER  IQ " » 19 »  ■  .  '  ) 

400  FORMAT<  1111  yTtO.  "THIS  IS  A  BLAST  PROBLEM”) 

600  PH R M. A  1(1111.1  4 0  ?  "  N 0  S U Oi  l  P R 0 BL E M  N U M B E R  -■  “  .13) 

700  FORMAT ( 124  ? "RADIATION  EFFECTS  HAVE  BEEN  INCORPORATED  AT  EARLY  T 
IS"  ) 

RETURN 

END 


SUBROUT I NE  CONTROL ( FLAG . TABLE » TSTOP ) 

DIMENSION  X  (  I  0 1 )  »U(  101 )  yPHOi  )  y  D  (  1. 0 1  ),Q(  101. )  y  FLAG  <  5  ) 
1  y TABLE ( 10) rCK(  7) y  DEBT ( 101 ) yEI ( 101 ) y  TEMP (101 ) yER< 101 > 


00 


I 

i 

i 

i 

i 


5  SO 
[,60 
570 


1  »  I  A  L I  F"  C  I  0  , 
n  im  iis  i: (i,! 

l.i  1.  i  1 1  NS  I.  f  t M 


v  U\  ( / )  v  i.i  i. :  i  1 1  <  i.  o  ,i.  j:  c  i  o.i  ),im  p  ( j.  o  j.  > ,  e  r  <  i  o  i  ) 
r  ( 2  y  J  o  i. ?  i. o ) .  in:: n  r  u-  <.  j. o i.  > ,  amass  <  j. o i  > 

1  0  1  )  v  H  I  AO  i  o  I.  ) 


'  ■:  :  ■  :  '  1  0  ’■  v , .'  RG:  LAP  <  I  0  )  >  <.■;'(  i  J  01  )  >  VC  <  iv 

'  *  '•  ‘  1  :  I  *  l  » 1 .  1  v  I  !  I  1 1  r  «  i :  . 

6  00  =  C - -  COLL  TO  OF  I  INITIAL  CONDI II Oil!!  - - - - 

610-  CALL  BLAST  (CAy  I  ALPHA?  Ml  y  LX  y  Rl  10  y  P  ?  D  ?  X  ?  0  L!  >  DT  ?  CK  ?  V  y  DELV ) 

61? 0  -  CALI.  RADSL  1  (  MAXNU  y  DT  ?  IJP T  fill'  y  TEMP  y  CK  (  2  )  y  D  ?  El  y  CK  ( 1  ) 

630-  1 «  LR  y  AMASS  -  X ) 

6  4  0  -  C - LI  S'  T  I N I T I A I...  CONDITIONS  - - - 

650-  TIME-CK ( 2 ) 

660-  IF  3 -FLAG (3) 

6  7  o  -  n  a  i.  u  I  c;  y  c:  i..  i ::  - 1.  ?  i  f  3 

6  SO--  A  A  -  M  OD  ( ICYCl  E  y  5  0 ) 

690-  IF  ( I  CYCLE  .  G  T  .  0  .  AND .  I.  CYCLE  .  LE  .  5  )  A  A =5 

700-  IF  <  AA . EG . 5 )  CALL  WRITER ( TEMP y  P  ?  ID  El t X  r DEBT  r ER ?  F ) 

710-  TFRAC-O.i 


7 2 0 -  CALL.  PLAN C K  <  N A X N U,  T E MP»B»W P T E M P ?  A V G N U ) 

730-  DO  A  NU=1 y MAXNU 

740=  UNITY-0. 

750-  DO  3  1=1 y 100 

760=  XXX-  FT  10 


:  770=  I F < I . EG . 1 )  XXX =0 . 3386 

j  780=  IF < UNITY. GE. 1 . )  GO  TO  3 

1  790=  UNITY-UNITY  t AMU <D< I ) y TEMP ( I ) ? 1 . ? XXX ) * ( X ( I + 1 > -X ( I  )  ) 

^  800=  IF<UNITY.GE, 1 . .OR. I .EG, 100)  JJ-IT 1 

i  810=  3  CONTINUE 

j  820=  4  RSFEAR(NU)  =  X(JJ) 

•  830=  DT-DTJ-DTH-1 ,E-6 

840=  I F  <  I C  Y  C  I..  E  .  L  T  .  1 0  0  )  D  T  =  1 .  E  -  7 

i  850—  IFCICYCLE.LE. 20)  DT-l.E-9 

;  860=  I F  <  FLAG <  5  >  . EG .  1  )  GO  TO  2 

,  870=  DO  1  1=1? 100 

j  880=  IF ( DEDTR ( I ) , EG . 0 )  GO  TO  1 

1  890=  I F ( D E LOCI). E G , 0 )  G 0  1 0  1 

|  900=  DT J=ADS <  TF R ACTE I ( I ) / <  DEDTR ( I ) *DELV ( I )  /AMASS  (I )  )  ) 

(  910=  IF(DTJ.LT.DT)  D  I'-DT  J 

j  920=  1  CONTINUE 

|  930=  2  CONTINUE 

940=  CALL  RAD FLO ( MAXNU ? X  ? DT ? DEDTR  y  TEMP ?  RSFEAR ?  D  ?  ER  ?  F  ?  D  y AVGNU ?  Rl  10 ) 

!  950=  DTH-'  DT  J=  l .  0 

i  960 -  I F  ( F L  A G  ( 5 )  .  I:  G  .  0 )  G 0  T 0  6 

I  970=  UG ( 1 ) =5 , 73 . 

{  980=  DO  5  I- 2?  1.00 

j  990=  5  VGCr. )  =  1.4 

i  1000=  6  CONTINUE 

i  1010=  DO  7  .1=1  ?  I  00 

«  1 020  =  A  =  D  (..))/(  VG  ( ,J  )  DP  ( ,J  )  ) 

.  1030  -  IF(A.LE.C)  GO  TO  7 

I  1040  -  D  T  J  ••  <  X  <  J  {  1  )  -  X  (  J  )  )  T  T  F  R  A  C  >)'  S  G  R  T  (  A  ) 

•'  1050=  IFCDT  J.LT.DTII)  DILI  D  I  J 

:  1060  7  CONriNUE 


1070  8  IF  (DTII.GT.DT)  Dill  HI 

1080  N  I  N  1  (  DT/DTI 1 1  I  1 

1090  Dill  DI/ILGAiCN) 

1100  DO  10  11  1 ?  i  I 


t 


01 


.  1100 
1110  c 
n : 

i 

I  1 4  0 

I I  5  0 

'  1 160' 
1170 
i  1180-- 
;  1190= 

j  1200 

S  1210= 

1  1220= 
i  1230= 

|  1240  = 
j  1.250  C 
|  1260= 
j  1270  = 
1280  = 
1290- C 
j  1300= 

|  131.0  = 

|  1320 =C 
!  1330= 

!  1340  =  C 
1  1350= 

•  1360= 

'  1370= 

1380  = 

»  1390= 
i  1400=8 
i  1410= 

|  1420  = 

(  1430= 

{  1440= 

!  1450= 

(  1460- C 
i  1470= 
j  1480= 
i  1490= 

|  1500= 

I  1510= 

|  1520  ' 
i  1530= 

!  1540= C 
j  1550= 

•  15A0=C 
i  1570  = 

!  1580= 

!  1590= 

:  i6oo 
!  1610 
1620  c: 

1.630  = 
1640 
I  1650= 


!UJ  10  II  1,N 

i..m  i  to  on'-1,  on  nor:  dmi  uii  p  - - - - - - - - 

1  .  .  =  :•  -  .  i!  '.i  .ii.r,/  !  i',1  i ‘i  i,'.  •  r:;i  in  y  Oi'T.i’iiT  *  m  i » iiT!  i  •  t  iii;-  y  v  i 

.  :  i  .  i  i  i  i  .1  •  v  I  i  ■: « !  v  i  >i  .i !.  '  •  y  ON  y  VU  y  V  i  i'!..I..O  i  08  ! 

.id  io  i  me.- 

m.  8 1  ci)  1 1  (i  1 1;«  ' . ; .;  < oik  :i:  >-  oo  ( i ) ) 
l  :i  <  i  i  i:.:t  <  j  )  ;  h  i  m  h  i  <  i  ) 
ir(Od).L'I.O)  11(1;  3.308C::  6 

tump  <  i: )  =tempii:r  <  ei  cm  temp ( i ) »  on <  i  >  > 

08(1 )  - 5 . / 3  . 

XXX  RIIO 

j|  (  I  .10.  I.  )  XXX  0.3386 
IS  1  ABUCi  1;  (  I  ) /l  .l  cm 

I  Fr  ( 1  5  I  .1  I  .  .01  )  08  v  I )  -  OCAMMA  (  EI  (  I )  y  D  ( I )  yXXX) 

10  P(  I  )  (  0  G  (  n-  I  >  1 1.  J  (  1  )  S1K  1  ) 

- CAI  l.  TO  81':  1  BOUNDARY  CONDITIONS  AND  TRACK  SHOCK - 

CAL  L  E.XI  AND  (HD  I  ALT 'll  A  *  N  i.  y  I.i  X  i  K'HO  1 1-' y  D » X  t  U  y  U  i  DT  y  CK y FL  AG  ( 5 )  ) 

IT  Ml:  CK  (  2  ) 

FLAC ( 5 ) =0 

. .  CALL  TO  COPY  PE  OOP  SI  ED  TABLES  - - - - 

i  f  ( r  j:  m f  .  c r.  .  i  c  i  o p .  o r  .  r  :i: me  .  o i: .  r a d i ... i :  ( i ) ) 

1 C  A I  l„  8  0  P  Y  (  X  >  IJ  1 1-  t  D  .  0  V  F I...  t .  G  (  2  )  y  I  I  h  E  ■/  T  A  D  L„  E  y  C  K  y  T  E  M  P  ) 

-  ADOPT  IF  TIME  LIMIT  REACHED  - 

I F  (  T 1 M  E.OE.T  S  T  OP)  R  E  T  l.l  R  N 

- . -  CALL  FOR  PE  ZONE  IF  REQUIRED  - 

IF  (Ml  .LT.94)  GO  TO  1.5 

C  A  L.  L  R  r:  Z  0  N  E  <  X  y  u  V  P  v  D  »  0  t  D  X » 0  A  M  M  A  t  ICY  C  L  E  y  T I M  E  y  N 1.  y  P II 0  ) 

PRINT*  y "AFTER  RE ZONE" 

c:  A  L.  L  W  R I T  E  R  (  T  E  M  P  » P  »  D  y  E I  y  X  y  D  E  D  T  «  E  R  y  F ) 

15  CONTINUE 

- -  CALL  RESULTS  AND  ABORT  PROGRAM  IF  CYCLE  LIMIT  REACHED  - 

CAI..  I...  C  0  P  Y  ( X  y  I J  y  P  y  D  v  0  y  F  L.  A  G  (  2  )  y  T I  i  i  E »  T  A  B  L  E  y  C  K  y  T  E  M  F:' ) 

PRINT* y " YOU  REACHED  THE  CYCLE  LIMIT  OF  "y I CYCLE 
P R I N I  * y “  TIM E  IS  1 y  T I M E » '  SEC." 

RETURN 

END 


S U B R 0 U T I NE  H Y D R 0  (  D  y  P »  X  y  0 . 0 E L  *  D X  v  C A  ,  I... ,  R H 0  y  G A M MA » Ni  r  DT  y  TIME »  F I. A G  t  C 
1 MAXNU  y  E I  y  DEDT  y  WP TEMP  y  TEMP  y  HR  y  F  y  AMASS  ?  ON  y  00  y  0  ,  DELO  .*  00 ) 

DIMENSION  DC  10...  )  y  PC  10:1.  )>X(  101  )  yQ<  101  )  -  L)  <  101 )  y  FLAG  ( 5 ) 

1  r  OG  (  1 0 1 )  y  OB.  <  1 0  I. )  y  AO  GNU  (10)  y  E I  ( 1 0 1  )  y  DEDT  (101).  00  ( 10 1 )  y  ON  ( 1 
1 0 1  >  y  T  E M  F  ( 1 0 1 ) y  R S F  E A  R ( 1 0 ) y D ( 1 0  I  y 1 0 ) v D E D T  R ( 1 0 1 ) , E  R <101) 

DIME N 8 1 0 N  F  <2.1 0 .1.  *10)  y  A M A S S  (  1 0  I.  )  v  D E I. R  <  .1. 0 :!  ) 

D 1  ME  NS  J  0  N  0  ( .1  0 1. )  .  h  E  L  0  (  1. 0  I. ) 

-  II  =  8  /  8  K  t  *  3  y  P = D  /  8 !  I  *  Sc  2  y  T  E  M  P = K  E  L  OIN*  E I  =  E  R = E  R  G  S/CM**3  y  00  =  0  N =CM**3/G 

INTEGER  FLAG 

-  INITIALIZE  LOCAL  PARAMETERS  - 

0  E  L.  <  1  )  0 
P I  =  A  C  0  S  <  - 1  .  ) 

QHAX-  - .  01 

DO  40  1=1 r 100 

IP  SORT  (  08  (  I  )  J  - 1 '  <  D/IK  I  )  ) 

. . -  --  8  0  M  P  U  T  L  N  r  IJ  0  r  S  8  0 8  I  i  I E  8  - 

Q  (  I)  0 . 5  TP  (  1  )  DP  1.0  ([)/’( I;  1  .1  <  X  <  I  )  T*  21  X  ( If  J.  )  t*2  )  TP  ) 

IF  (  (  OI  L  <  1.  I  1  )  OI  L.  (  I  )  )  .  GL  .  0  .  )  Q  ( I )  =0  . 

40  CONTINUE 
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I. 650 
1660 
I  670 
I 

169  0 

J.  700  = 
17  10 

J  7  2  0  “  ■  C 
1730  = 
1740  = 

1 750  = 
1760  = 
1770  •(: 
1700- 
1790:= 
1800=: 
1810  = 
1820- 
1830= 
18-10  = 

1  830  = 
1860  = 
1870  = 
1880  = 
1890  = 
1900  = 
1910  = 
1920  = 
1930= 
1940  = 
1950  = 
1960  = 
1970= 
1980=  C 
1990  = 
2000  = 
2010  = 
2020-C 
2030  = 
2040  = 
2050  = 
2060  = 
2070= 
2080  = 
2090  = 
2100  = 
2.1.1  o  = 
2120  - 
2130=8 
2140  = 
2130  = 
2160  = 
2170  = 
•’180  = 

.  190- 
2200= 


4  0  ..utii  :i » aji : 
BEl.RC  l  )  0. 

00  4  2  1  1x99 


'  •  •  -  ■ !  1  '  •  I  ■  1  ’  -i  ■  J  M  1  ;  /  u  )|  I Ii :> ( i >  : ,  mA:>b  <  j,  t  j.  >  > 

V.  !.u  ;  i  >  ole ( j. ;  ,i  )  ; ,  .{.c;  in 

4 2  1  ■  8. L. ft  ( .1  {  :l  )  =  1.1 1  ;.  (  0 i. I  ( .1  f  J.  )  1  0 . 5 X  A C C ,{; I.i T  ) 

I'D  45  I  =  1x99 

. .  ASSUMES  6,.!  ft  I:  III!  A  Oft  8  AS  AN  1 DEAI  GA8  DURING  RAIil  AT  I  ON  i'll,*.:  f 

i::i  ( I  >  cl  <  I )  I  4  .  *PI.;XP(  I )  ru  ( I )  )  *  (X(  I )  *.?.2*DELR  <I)-X<I  +  1>  **2S.I'it"LR  (  T-i 
1 )  /  AN AGS ( I ) 

X(IM)  =  X ( 1 M.  >  I-  D E L ft ( 1  1 1  > 

45  CONTINUE 

- - LOCATE  NEW  SHOCK  POSITION  - 

DO  50  .1  =  1x100 


OCX  I )  =  1  ./IXJ ) 

oo  ld =0  <  :i:  > 


0(1)  =4  .  /3  .  >:<p I  XX  X  (  I  i  1  )  9 * 3 -  X  (  I )  9*3  ) 
IX  I )  =  AMASS  ( I )  /V  ( 1 ) 

DELO  ( I )  •:  ADS  (0(1)  -OGL.D  ) 

ON ( I ) =1 ./IX  I) 

50  CONTINUE 

DO  52  1=1x100 

IF ( UMAX. GT. 0(1) >00  TO  52 


QMAX--0  ( I ) 
N 1 = I +1 
52  CONTINUE 
RETURN 


END 

SUBROUT I NE  BEAST ( CA  x L  x  N 1 x  DX  x  El  10  x P  x  D  x  X  x  Q  x  U  x DT x  CK x  0  x  DELO ) 


DIMENSION  P(lOl)  fix  101)  xXCI.01.)  xQdOl)  xU(lOl)  ,CK(7)  xO(lOl)  »  DELO  <10.1 
1 ) 

G  A  =  4 

P I  =  AC  OS  (  -1. .  ) 

-  READ  YIELD « SCALING  PARAMETERS x AND  PROBLEM  GEOMETRY  - 

READ*  x  YK I x  SP x  SD  x S  T 


L.  =  3 

UR  I TE* ,  YKT  x  SP  x  SD  x  ST  v  L. 

-  COMPUTE  INITIAL  CONDITIONS  USING  SCALING  LAWS  - 

P N 0 R M  - 1 .01 4E6*SP 
ck(2)=t :mr:= .1.  ,e-7 
PBLAST=9.8E 13 
DX  8V. 

C  K ( 1 )  =  YKT 
RIIG  1 . 293E-  3TSD 
N 1  =  1 

CK  (  5  )  =  P 'NORM 

P  R  I  N  r  X;  x  "  I '  N  0  R  H  =  *  x  P  N  0 1  ( M  x  *  P  D  L  A  G  I  =  ’  x  P  D  L  A  S  T 

PRINT*,  “  T  I  NE  = "  x  TIME  x  "  DX  =  *  x  D>:  x  “  Rl  IIP  "  x  Rl  10  x  "  Nl  =  *xNl 

-  AS SION  SCALED  VALUES  TO  PROBLEM  ARRAYS  - - - 

DC)  1.0  1  =  1x101 


CX  I  )  ••  T  i  (  1  )  0 
IK  I  ) =RHO 

if  ( i .  eq  .  i ) x  ( :i: )  :  0 

IE  (I  .or.  I.  )X(  I  )  X  (  I  I  )  I  DX 
1F(  J  .1  E  .  N.I  )P  (  I  )  P)  I  AST 
IF(  I  .  0  T  .  N 1  )P(  i  )  -I'NORM 


X — 
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2? 00  -  I F ( J . 0 T . N 1 ) P ( 1 )  - P N 0 R M 

:>:>:!  o  ■  i  o  re  i  i,'  (  t  >  -  o . 

’  >  !'  |  ;■  l  ■ 

■'  V  (  L  )  ■  ♦  /  x:>  <  I  ■*  1  '■'!■■  ( /\  (  .1.  T  .L  ) 3  *"■  X  v  J.  )  $  2* ) 

2',Mo  20  con  (  inwl 

2250”  CK<3>  «  CPUS  *  X<2) 

2260  -  RET  URN 

2270-  END 

2200-C 

2 2 9 0  -  B U D R G l)  I  I N E  E X P A N D  (  CA  y  L  y  N 1.  vPXy  R I  I 0  y  P  i  D  »  X  •  Q  y  IJ  y  D T  y  C K  y  F I... A G  ) 

2300-  It  I  MENU  I  ON  P  <  1 0 1  )  y  D  ( 1 0 1 )  y  X  (  1. 0 .1.  )  y  G  ( 1 0 1 )  ,U(  101)  r  CK  (  7  ) 

2310-  INTEGER  FLAG 

2  3  2 0  -  I F  ( F  L  A  G  .  EG  .  1.  )  II  M  I T -  K I  IF  TV  -  C  Y  C  L  E  -  0 

2330  C -  RESET  SYMETE.IC  B.  C.  AND  ADVANCE  TIME  - 

2340-  X(l)-U(l)-0. 

2  3  5  0  «  I F'  ( F  L  A  G  .  F.  G  .  I.  )  P  N  0  R  M  -  C  K  (  5  ) 

2  3  6  0  •  =  I F  <  fr  I...  A  G  *  E  0 . 1 )  T I M  E  -  C  K  ( 2 ) 

2370 -C -  EGG ATE  NEW  SHOCK  POSITION  - 

2300-  IP-N1 

2370 -C- TIME  SHOCK  PROGRESS  FOR  SO  ITERATIONS  - 

2400~  KFII  TY-K FIFTY  I  1 

2410-  TIMI1 IIMlT  i  DT 

2420-C - -  COMPUTE  OVERPRESSURE  AND  OVER DENSITY  - 

2430-  PM  A  X  -=  D  M  A  X  -  V  h  A  X  •  0  . 

2440-  DO  20  I-J.  y.1.00 

2450-  I F  ( l,  ( I )  .  G  T  .  V  M  A  X )  V  H  A  X  -  IJ  ( I ) 

2460-  I F  ( P  <  I )  .  G  T  .  P  M  A  X  )  P  M  A  X  -  P  ( I ) 

2  4  7  0  -  I F  <  D  (  I  )  ,  0  I  .  P  h  A  X )  h  ,s  i  A  X  -  D  <  1 ) 

2480-  20  CONTINUE 

2  4  V  0  ==  C  K  ( S  )  =•:  OP  -  p  M  A  X  -  P  N  0  R  M 

2 500-  CK  < 6 ) - 0 D - D M A X - R 1 1 0 

2510-C . . —  PLACE  COMPUTED  VALUES  IN  CHECK  (CK)  ARRAY  - 

2520*  CK  <  1 )  -CYCLE -C  YCLL-I  1  . 

2530-  CK ( 2 ) -  I  1  ML-  I  I  MET DT 

2540-  CK(3)-DT 

2 5 5 0 -  CK  <  4  )  * P 0 S X  ( 1 P ) 

2560-C . —  COMP HIE  NEW  SHOCK  MATERIAL  VELOCITY  - 

2570=  CK ( 7 ) =  V  M  A  X 

2580-  IF (KF1FTY.LT. 50)  GO  TO  30 

2590-  OPUS-FOB 

2600=  TIMIT-K FIFTY-0 

2610-C - -  WRITE  CHECK  VALUES  EVERY  50  ITERATIONS  - 

2 6  2  0  ---  U  R I T  E  (By  2  0  0  0 )  (  C  K  ( I )  .  I  -  J.  y  7  ) 

2630=  IJR  II  E  (  2  -  :l  00  )  (  CK  (KG  K-l  y  7  ) 

2 6  4  0 -  1 0 0  I  ■  0  R  M A T ( 2 X  y “ C Y  C I  E  " y F 6 . 0  t  *  TIM E - " y 1 P E 9 . 3  y  *  D T - “ • F 9 . 3  y 
2650-  1*  S  I 'OS  -  * y  E9 . 3  y "  0  PRES- " y E9 . 3 y  “  0  DEN- * E7 . 3 y “  S  VEL=* 

2660-  2  y  E  V . 3 ) 


26/0  = 
2680 
2670- 
2/00 


2000  F ORMAT  <  I  X , n  CYCLE  - " y T  5 . 0 y ■  T ( SEC ) - ■ 

.1  y 7 . o y  *  i tis (CM)  " . i: v .:v  o  1  ■  ( u / c m * t 

2 ft::?. 2 y  ••  m v i I.  < l; m / s r: c ) ---•  •  -ev  . 2 ) 

30  RT TURN 


1PE9.2» ’  D ! (SEC)-" 

>■-  " »  e 9  •  2.  y  “  n:  a;/ cm« 3) 


2710  END 

2720- C 

2730  SIIPK-OI I  I  1  Nl .  RF  70NI  .  (  X  .  U  y  P  y  D  ,  0  .  DX  .  GAMMA  v  1  CYC!  I  v  I  IMF  *  N I  »  Rl  10  ) 

2740-  D  I  Ml  NS  I  UN  X(  101  )  ,U(  L01  >  ,PU  01  >  .  D(  I  01  >  yPUD(  .1  01  y  3  )  rO(  10 1  ) 

2750=  Nl-NL/2 
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950 

,o... 

N  !  = 

1  };•  T 

=  N  .1 

r  r  ; 

•'  70 

V.1  t  \  \ 

Du 

L  1  1 

1 0 

2/uo= 
2V9o-C 
2800 -■■■■■ 
28.1  0 
2020== 

2  0  3  0  == 
2840-C- 
2830- 
2060  =  : 
2870- 
2880- 
2890-C- 
2900 -- 
2910" 
2920-C- 
2930" 
2940== 
2950= 
2960  = 
2970- 
2900- 
2990==C- 
3000- 
30.1 0  =  0- 
3020- 
3030- 
3040-= 
3050  - 
3060- C 
3070== 
3000-= 
3090-C- 
3100== 
3110== 
3120--C 
3130== 
3140  = 
3150  - 
3160== 
3170- 

3ioo=n- 

3190- 
3200  = 
3210  C 
3220== 
3230- 
3240== 
3250-= 
3260 
3270  = 
3200=  c: 
3290  = 
3300== 


2  *  1  00  )  T  I'M  v  ? :  Y(  I  . 

J-  !.  y  99.2 
1=  <81  1  )/2 

EIJ1==EIA(P(  J>  y  1 1  (  J  )  ;  R 1 1 0  ) 

E  I J2-E I A  (  P  (  J  !  1 )  y  1.1  (  J  i  1 )  y  EMO  ) 

A G A ii A  =•  0 G  A I  i  A  (  L I J  j  ,  0  <  J  )  R |-| 0  ) 

DGAMA==VGAMMA  ( P  I J2  y  D  <  J  !  1 )  y  I •; 1 1 ( j ) 

. .  COMPUTE  MAGS  AMD  NEW  DENSITY  - 

XMJ 1  =  D  (  J )  >!•  (  X  ( ,H  1  )  X*3-  X  (  J  )  «3  ) 

XMJ2-D  ( J }  1 )  ;i:  (  X'  (  J  f  2 )  «3  -X  (  J  I  1  )  XX3 ) 

PUD  ( I  y  3 )  ==  (  XM  Jl  -i  XMJ 2 )  /  ( X  (  J  !  2 )  **3-  X  (  J )  «3  ) 

X  h  1 1  ■-=  P  U  D  <  I  y  J  )  :l:  (  X  (  J  i  2  )  *  *  3  -  X  (  J  )  T 1 3 ) 

-  COMPUTE  OLD  CELL  CENTERED  VELOCITY  - 

U  Jl  =  ( 0 ( J ) -MX  J!  1  )  )/2. 

U  J2-  <  IJ  ( J 11 )  {  U  (  J  12 )  )  / 2  . 

-  COMPUTE  OLD  ENERGIES  (.  KINETIC »  INTERNAL  y  TOTAL) - 

EKJ1  ==UJl**2/2. 

EKJ2==UJ2**2/2. 

E I  Jl  -P  <  J )  /  <  (  AG  AM  A- 1 .  )  *D  ( J )  ) 

El  J2-P  (  J+l )  /  (  (  DO  AM  A  - 1 .  )  >XD  (  J  M )  ) 

ETJ1-EKJ1  i-EIJl 
ET  J2--EK  J2F  !I  I J2 

-  COMPUTE  NEW  TOTAL  ENERGY  - 

E  T 1 .1.  (  X  M  J 1  *  E  T  .J  .1.  X  M  J  2 :}:  E  T  J  2  )  /  X  M 1 1 

-  COMPUTE  MOMENTUM  - - - 

XMOMJ 1.  =-XM  J 1*UJ 1 

XM0MJ2=XMJ2*l.JJ2 

PUD  <  I  y  2  )  --=  (  XMOMJ  1  fXMOM J2 )  /X M 1 1 

XM0MI1.==XMI  1  *PUD<  I  »2) 

-  COMPUTE  NEW  KINETIC  AND  INTERNAL  ENERGY  - 

EKIl-PUIK  I  »2)**2/2. 

El  1 1  l-ETI  1  ~EI\I  1 

GAMMA-VGAMMA  <  E 1 1 1  y  PUD  ( I  y  3 )  v  RI  IO  ) 

P  U  D  ( I  y  1  )  =  E 1 1  J.  J  (  G  A  M  M  A  -  J  .  )  >;<  P  U  D  (  I  »  3  > 

---  BYPASS  LENGTHY  REZONE  PRINTOUT  - 

IF<N1 .GT. 10) GO  TO  10 

WR I  T  E  (  2  »  200  )  X  (  J  >  r  U  Jl  i  P  (  J  )  y  D  (  J  )  y  XMJ1  yXMOMJ.1.  y  ETJ1  y  EK  Jl  y  E  I  J  1 


(J!l  )  *  U.J2  y  P  (  J  {•  1  )  y 
<  J)  yiJJJ  V  PUD  <  J  y  1  )  y 


D  (  '  i 
PUD  ( 


VELOCITY  FOR  CENTER  AND  RESET  DX 


VALUES  IN  PROGRAM  ARRAYS  - 


WRITE  (2y  200): 

WRITE ( 2  y  250 ) : 

10  CONTINUE 

-  ASSIGN 

PUD  (  1.  y  2 )  0 

DX==2  .  T:  D  X 

-  STORE  NEW 

DO  20  1=1 r 49 
J  2 :!  I  - 1 
X ( 1 ) -X ( J ) 

P ( I )  PUD ( I y 1 ) 

II<  I)  PUD  (1.3) 

IF  (  I  . PO . 1 ) OUT  0  20 

-  COMPUTE  NEW  BOUNDARY  VELOCITIES 

0  (  I  )  =  (  PUD  <  I  I  y  2  )  t  PUD  (  .1  y  7  )  )/2  . 

20  CONTINUE 


I )  y  X  M  J  2  »  X  M  0 !  1 J  2  v  E  T  J  2  » CL  K  J  2  y  E I 
3  )  i  XM 1 1 »  Xi  tOM  1 1  y  ET  1 1  y  EK  1 1  » 
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3300  r 
3  V  '  ■ 

.*,330 

.  3340 
3  3 ,  >  {/  ■  ■ 
3360  = 
3370= 
3380=C- 
3390= 

;  3400= 

:  3410= 

:  3420  = 

\  3430=C- 
;  3440= 

;  3450  = 

!  3460= 
i  3470= 

;  3400= 

:  3490= 
i  3500= 

>:  3510  = 

|  3520= 

■  3530= 

;  3540= 

3550  = 

■  3560  - 
3570  = 

■  3500= 

}  3590= 

(  3600=5 
:  3610  = 
3620  = 

I  3630= 
j  3640= 

/  3650- 
j  3660= 

’  3670= 
3600  = 
j  3690= 

!  3700= 

,  3710=0- 
j  3720= 
i  3730= 

|  3740= 

J  3750= 

!  3760= 

1  3770= 

:  3700  = 

■  3790= 

'  3000 

3010  = 

3020 

3030 

!  3040  0- 
:  3050= 


APR AYS  - 


20  CONTINUE 

- .  on- 1 mo:  051111:  0  j  or  upper  half  of  arrays  - 

D  LI  .  j  .1.  1  i  \f  >■  1  1  /  .i. 

X  <  1.  )  X  (  4  -1  MUX 
li(  1)0(1  )  0 
PCI  )  O'- ( 101  ) 

1.1  <  1 )  ncioi  > 

25  continue 

-  COMPUTE  NEl-J  VISCOSITY  - 

DO  30  1=1 y 50 

0(1)  =0  .  *I.i  <  I )  *  <  U  (I  T  1 )  - 1.)  ( I )  )  *  *  2 
I F  (  U  <  I T 1  )  -  U  ( I  )  .  Cl.. .  0  )  Q  (  I )  =0 
30  CONTINUE 

-  SET  HOUND ARY  CONDITIONS  - 

U< 1 )=X< 1 ) =0 
F ( 1 ) =P ( 2 ) 

0(1) =0 ( 2 ) 

D  < 1 ) =D ( 2 ) 

i.  0  0  F  0  R  H  A  T  <  1.  II 1  y  T  2  0  y  ”  R  E  7 0  N  E  T I M  E  =  "  ?  G  9 . 3  y  "  IT  E  R  AT  I  ON  N IJ  M  B  E  R  ’  ,  1 7  , 
1  /  /  /  /  y  T  4  y  " P  0  s  I T 1 0  N  "  y  T  :l.  iif  "V  E  I...  0  C I T  Y  *  »  T  3  2  y  "  P  R  E  S  S  U  R  E  *  y  T  4  6  y  *  D  E  N  SIT  Y 
2T62» *  MASS " »  T74 1  " M 0 M E NT U H " y  T 8 8 y " I  E N E R G Y " y T 1 0 2 y * K  E N E R 0 Y " y 
3  T 1 1 6  y  "  I  E  N  E  R  G  Y'3/,1 7  y  “  C  M  "  y T 1 9  y  “  C  M  /  0  E  C "  y  T  3  3  y  *  D  /  C  N  *  *  3  ■  y 
4 T 4 6  y  ■  G/CM  Z  "  y  T63  y  "  GM  “  y  T76  y  0  G  -SEC  "  y  T9 1  y  *  ERG/G  "  y  T 1 05  y  ■  ERG/G  •  y 
5T119- “ERG/G" y ///) 

2  0  0  F  0  R  M  A  T  ( 3  X  y  9  (  G  9 , 3 »  5  X  )  ) 

1 5  0  F  0  R  M  A  T  <  5  X  -  9  <  G  9 . 3  y  5  X  )  y  /  ) 

RETURN 

END 

r  U  N  C  T 1 0  N  V  G  A  M  M  A  (  E  I»Dr  R II 0  ) 

- -  HUFF  ROUGH  VARIABLE  GAMMA  COMPUTATION  - 

E1E-EI/1.E10 

IF (El E. GT. . 131 5) GO  TO  10 
GMONE-- .  3981 
GO  TO  30 
10  CONTINUE 

IF  (  EIE  •  GT  .  1. .  )  GO  TO  20 
G M  0  N E  =  - . 0 3  9  9 * (EIE - ,1315) * *  2 ! . 3  9  3 1 
DO  TO  30 
20  CONTINUE 

G M 0 N E  .  .1 6  +  ,62  4/(2.  I  E I E  ) 

-  adjustment  for  density  - 

RIIOX  =  D/RIIO 
EL=AI.OG  10  (  L  IE) 

AL  PI  IA-- .  0 f 1 7 7 T F I...  .  OPlOt  El  **2  +  .  0035*El„**3 1  .  0002:!  :"i.  T * 4 
GNOME  ■  OMONI  •:  IN  IOZTTAI  .F  I IA 
30  VOAihiA  OiiUNP  !  1  . 

RETURN 
t.  ND 


siiDRoin  ini:  copy ( x >•  id p , d , o, n  ad « time . tai  i  >  l. 

D I  Ml  NS  ION  X  i  1  0  I.  )  y  I.J  (  I  0  J  )  ?  1 '  ( 1  C  L  )  y  D  l  i  0  I  )  y  u  <  |  0  1 
1  y  CK  (  7  )  y  II  1 1 1  ’  (  i  01  ) 
imiTTi;  n  ao 

-  RF-.SU  CAl  LING  PARAMF  I  FKS  - . . . 

F  l.  AG-f  LAG  -  .1 


.  .10) 


1 


3! ,n  - 


C 


3KNV 
3S:">' 
o  vc  o  • 

391  o 

3970= 

3930- 

3940= 

3VS(>-; 

3900 

3970  = 

3900= 

3990  = 

4000  = 

4010= 

4020= 

40  30- 
4040  = 
4000  = 
4040  = 
4070  = 
4080- 
4090  = 
4100  = 
4110= 
4120  = 
4130  = 
4140=- 
4U>0= 
4190  = 
4170  = 

4 1  80  = 
4190  = 
4200  = 
4210  = 
4220  = 
4230  = 
4240  = 
4280  = 
4240  = 
4270= 
4200  = 
4290  = 
4300  = 
4310- 
4320  = 
4330  - 
4340  = 
432)0-- 
4340  = 
43  70  = 
4300  = 
4390- 
4400  = 


RAG  I  I  00  I 

Hi)  )(.  I  y 

I  v  I  Oil :l  1 .  <  I  )  1  ADEL  USD 

t,r  (( i  I'tit . i  i  .oraki  i  <  i  >= 1000. 

- -  WRITE  INFORMATION  ON  LOCAL  FILE  (FILES)  - 

CNii  v  i.i. ory 

worn  < 7 ...  i  coomnc,  <  x  <  :i  >  y  u  <  x )  ,pa  ))D(I).h(d»  i=i  rod 

URI 7  r:  (  7 . 2000 )  (  ok  <  i; )  ,  Iv  1 ,7) 

WRITE  ( 2  ,  1.0)  0 )  T  I  hr 
CALL  I 'LIU  SM(X,P,  100) 

U R 1 1  [:  <  2 >  2 000)  (CLCK)  »  K  =  1,7) 

WRITE  <2,  i040)TlME 

call  i  lot  so  < x  ,  temp  » .100 ) 
ur  n  r  ( 2 , 7000 )  <  ck  (  k  >  ,  i<~ if?) 

WRITE (2» J 020) TIME 
CALL  PLmSfKXf  rnlOO) 

UR I TE ( 2  *  2000 ) ( CK  C K ) > K  =  1 , 7 ) 

UR  1 7  E  <  2  v  :l  001 )  TIME 
CALL  PLOTS!  i  <  X  ,  U  ,  1 00 ) 

UR  I TE <  3 ,  2000 ) <  CK <  K ) , K  = 1 ,  7 ) 

WRITE ( 2 f 1030) TIME 
CALL  PLO(EiKX,Q,  TOO) 

W  R I T  E  (  2  ,  2  0  0  0 )  < C.  K  ( K  )  ,  K  =  1,7) 

1000  FORMAT (IX, "TIME  IS" »F9.5, "SEC" ,//» T7» " POSITION  * ,T22 

1 V  •VELOCITY"  ,  T37  ,  "  l:'  RE  SB  LIRE  *  » IS  2  ,  “DENSITY"  ,T67»  "VISCOSITY" 

2 »  /  >  T 9 *  «  (CM)  *  ,722,  "  (  CM/ SEC)  "  ,T37,  "  <D/CM**2>  *  »T51 1  ’  <  G/CMTT3  )  * 
3,767,  "  ( li/CM-f  I  2  )”»//»  10:1.  ( ’3  C  0  X  ,  IFEU),3  ?,/>  > 

1001  FORMAT  ( .1111  rT4()»  “VELOCITY  (CM/SEC)  VC  POSITION  (CM)  AT  TIME  =  ' 
1 » I PE 10 *3- *  SEC" ) 

1010  FORMAT  (  1  i-l  1  .  TOO  *  "  PRESSURE  ( D/CMTY2.  )  VS  POSITION  ( CM )  AT  TIME  = 

1 f IP El 0.3, *  SEC " ) 

1 0  2 0  F 0 R M A  T ( 1 1 1 1 »T 4 0 , " D E N 3 1 T  Y(G/CM * 1 3  )  VS  POSITION(CM)  AT  TIME  =  ’ 
1 lPE10.3v  *  SEC" ) 

1030  FORMAT  ( 1.111 »  T40»  "  VISCOSITY  ( H/CM**2)  VS  POSITION  (CM)  AT  TIME  - 
1 v 1 PE 10. 3v ”  SEC “ ) 

1040  FORMAT (1U1 »  T 40  » " TEMPERATURE < KELVIN )  VS  POSITION (CM)  AT  TIME 
1 V  :l.  PE. 1.0.3*  “  SEC"  ) 

2000  FORMAT  ( IX,  "  CYCLER-' "  *F5.0,  "  T<SEC)  =  *  ,  1PE9.2*  *  DT  (SEC>«* 

1  ,  h. 7 . 0  *  ”  POS <  CM )  =  "  *  E 9 . 2  *  "  OP  ( D/CMTME  )  =  *  ,  E9 . 2 ,  "  OD  <  G/CM**3 )  “  ’ 
2 , E 9 . 2  v  *  MVEL C CM/SEC )  =  *  *  E9 . 2 ) 

RETURN 

END 

SUDROUT I  Ml  PLOTSM ( X  v  Y  v  N ) 

DIMENSION  X<100)  ,Y(  100)  ,W(i(>0»40)  v  WL  (  40  > 

Nf'l  ■-  N  T  I 

DATA  SLASH  *  DOT,  STAR,  OH  v  DRANK/ l.H/v  1H.  r  II 'Tv  1H0*  "  *  / 

XflAX- X  (  N  I  J  ) 

YMAX  YM] N  0 
DO  10  I  1  v  N 

IRYd).  OT  .  YMAX  1 YMAX  Y  ( I ) 

IT  <Y(  I)  .1  I  .  Y  M  t  N  )  Y  M  I  N  Y<  I) 

1.0  CNN  l  INI  h: 

YOVt  K  YMAX- Y ( 1 00 ) 
no  30  1  L  V  : 1.00 
DO  20  J  1,40 
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44<)0~  DU  20  J—  1*40 

4410-=  IK  I.  J)= BLANK 

4420"-  •  r  <  t  .  rn .  i.  hk  i  » j>  si  ash 

4  4.-,')  ',1.1  >  i  J  ■,  I  ■■  J  ■  VO  V 

44  40  -  W (10* 20 ;  U ( 20  *  20 ,  W < JO » 20 ) =  W (40* 20 ) ( SO  *  20 ) ~U ( 99 » 20 )  =  *  0 ' 

4 4 S 0 U  ( 6 0  *  2 0 )  W  (70*2 0  >  •  W  (  8 0  *20)  =  W  (  90  *  20 )  =  W  ( 1 00  *  20 )  =  *  0  * 

4460  =  i.J  ( 9  *  20  )  U  (  98  *  20 )  "  1  “ 

4470-  W  ( 1 9  *  20 ) -  *  2 “ 

4400  •  W  (29  *20)==*  3* 

:  4490-  W  <  3  9  *  20  > -  *  4 " 

■  4500  =  W ( 49  *  20 ) - " S " 

■  4510-  W  <  59  *  20 ) -  *  6 ” 

i  4520=  U (69 *20)=" 7* 

|  4530=  W  ( 79  *  20 )  = " 8  * 

:  4540=  W ( 89  *  20 )  =  " 9 " 

1  4550=  20  CONTINUE 

I  4560=  SCAl.E=YMAX/20 . 

f  4570=  S T 0 R E  =  - Y M I N / 2 0 ♦ 

4580=  IF  ( SCALE  ♦  LT  .  STORE )  SCALE-STORE 

4590=  DX-XMAX/5 . 

1  4600=  SCA!J>SCALEi.  UPSCALE 

;  4610=  G M A X = 2 0 ♦ £ S 0 A L E 

4620=  DO  30  1=1*21 

?  4630=  X I M 1  -  ( I  - 1 ) 

;  4640=  8 T 0 R E = X I M 1 1 SCALE 

•  4650=  K=21-I 

I  4660=  l..  =  19  f I 

4670=  WL  ( L  )  =- STORE 

4680=  WL.  ( K )  -STORE 

4690=  30  CONTINUE 

!  4700=  DO  40  1=1*26*5 

4710=  J=4*< 1-1 ) 

:  4720=  XI=( 1-1 )/5 

:  4730=  UL<I>=XI*DX 

4740=  40  CONTINUE 

i  4750=  D X X  =  X M A X /  .1 0 0  . 

1  4 760=  DO  50  1=1 *N 

j  4770=  DO  44  K=1*N 

j  4780=  XK=K 

!  4790=  IF  ( X  ( I. )  .  GT  .  XK#DXX-DXX/2 .  .  AND  .  X  ( I )  .  L.E  .  XK*BXX+DXX/2  .  )  GO  TO  45 

i  4800=  44  CONTINUE 

i  4810=  45  CONTINUE 

j  4820=  DO  50  J  = 1*40 

1  4830=  XJM1= J-1 

’  4840=  XJ=J 

!  4050=  IF  (  ( GMAX-X  JM 1 T' SCALE  . GT . Y ( I )  )  .AND.  (GMAX-X J#SCALE 

4060=  1 .  LE  .  Y  ( I )  )  )  W  ( K  *  J )  =OTI 

;  4870=  50  CONTINUE 

:  4080=  WRITE ( 2  *  J  00 ) XMAX *  YHAX *  YMIN  * YOOER 

\  4890=  K-5 

4900=  WR I TE ( 2  *  200 ) DM AX 

'  4910=  DO  60  J  -1 .40 

4920=  If  ( J.EiJ.KHU)  TO  55 

4930  =  WRITE  (  2  »■  TOO )  (  W  ( I  *  J )  *  1  =  1  *  N  > 

4940=  01)  TO  6  0 

:  4950=  55  K- K \ 5 


’Tr" 
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4950  =  55  E-Ki5 

4  9  6  0  ~  ( 1  !■'  I  T  I-  <  2  /  2  0  0 )  W  L  (  J  )  z  (  W  (JtJ)t  J  =  j  ,  M  ) 

■v:  ■■ 

^ i'--  .•  ;  .n i.'cs  zc-o)  <wl < i )  > :i  -  i,2o/() 

4 f/ V o :■  1 00  I  t) R MAT  ( ///  /  T 3 0  /  "  1 1 1 J  S  p L D T  I ! A S  X M A X -  "  ,  I.  p El  0.3/  *  >  YMAX-  * 

5 000-  1  - E 1. 0 . 3  »  “  /  V M  I  N •••■=  11  .  E  .1. 0 . 3  .•  “  A N H  Y  D V E R  -  “  ,  E 1 0 . 3 ,  ///  > 

5  o  1 0  -  2  0  0  r  G IV  r ,  A  T  ( 5  X  v  G 1 0.3/  2  X »  I  0  0  A 1  ) 

5020-  300  FORMAT  ( /// »  T  1  4  »  A  <  :l. PC  1 0.3/ 1  OX )  ) 

5030-  400  FOIVMATd/XzlOOAl) 

5040-  RETURN 

5050-  END 

5060-C 

5070=  BUD  ROUT  I  ME  RADI  LC)  (  MAX  MU z R z DT  /  DEDT  /  TEMP  /  RSFEAR  z  B  /  ER  /  F  /  D  /  AVGNU  /  RHO ) 

5080-  DIMENSION  R(iOl) / RSFEAR ( 10 ) / F ( 2 / 101 » 10 ) /DEDT(lOl) 

5 0  9 0 -  D I M  E  N  S 1 0  N  D  ( 1 0 1 )  /  A  0  G  N I J  (10) 

5 1 00-  D I  MENS  I  ON  B  ( 1 0 1.  » 1 0 )  z  TEMP  <  1 0 .1  >  /  ER  (101) 

51 10=C -  R-CM  /  TEMP-RELOIN ,  ER  -ERGS/CNTT3  /  DEDT -ERGS/  (  CM**3-SEC ) - 

5120-C -  F-B-ERGS/ ( CM«2-SEC )  /  RSFEAR- 1 /CM  - 

5130-  P I -AC 03 ( -  1  .  ) 

5140-  DO  30  NU  -  1/MAXNlJ 

5150-  F  ( 1  / 101  /  NL) )  -  0. 

51 AO-  F ( 1 / 1 / NU )  -  0. 

5170=  F ( 2 / 1 / NU )  -  0. 

5180-  DO  15  J  -  1/99 

5190-  I  -  101  -  J 

5200-  FIB-PI ( I / NU ) 

5210-  XXX -RHO  : 

5220-  RH1-  R  < I  + 1 )  -  R ( I ) 

5230-  TAU-RM  IT  AMU  (  D  ( I. )  /  TEMP  ( I )  ,  AVGNU  ( NU  >  /  RHO  ) 

5240=  UMAX-  TAU  *  SORT  (  ( R  ( I  1 1 ) -!  R  ( I. )  )  /  RM1 ) 

5250-  IF (I. EG. 2)  XXX-0.333A 

5260-  TAUM-  <  R  <  I )  -  R  (I  -1  )  )  *  AMU  ( D  ( I  -- 1 )  z  TEMP  ( I  - 1 )  » 1 .  /  XXX ) 

5270-C  I F  (  TAU .  GT .  2 .  .  AND  .  TAIJM  .  GT  .  2  .  )  GO  TO  5 

5280-  GO  TO  10 

5290-  5  Fd/I/NU)  =ABS(8.*PJ  *  (B(I-lzNU)  -  B(IzNU))  /(3.*<TAU  +  TAUM))) 

5300-  GO  TO  15 

5310-  10  Fd/I/NU)  -  PIB  f(FdzldzNU)  PIB)  *  <<<R(I  +  1)+R(I>)  *TAU 

5320-  1 ) TT2  /  <4.  *  R(I)**2)  *  (G(TAU)  -  G(UMAX))  -  RM1**2  /  (2.* 

5330-  2 ( R <  I. )  *  TAU)«2>  t  (II (TAU)  -•  ll(UMAX))) 

5340-  15  CONTINUE 

5350-  BO  30  I  -  1/100 

53A0-  XXX  —  F(  II 0 

5370-  RS2  -  RSFEAR (NU)  **  2 

5380-  RM1  -  R(I+1)  -  R ( I ) 

5390-  IF ( I . EQ . 1 )  XXX-0.338A 

5  *  0  0  -  T  A  U  -  R  M 1  *  A  it  U  (  D  ( I )  /  T  E  M  P  ( I )  z  A  V  G  N I J  ( N  U )  /  X  X  X ) 

5410-  IF  (I  .ME.  100)  T  A  UP  -  ( R  ( I  +2 )  --R  ( 1 1  1 )  )  TAMIJ  (  D(  HI )  /  TEMP  ( 1  +  1 )  / 1 .  /  RHO) 

S420--C  I F  <  TAU  .  A  t .  2  ,  .  AMD  .  TAUP .  GT .  2  .  .  AND  .  I .  NE  .  1 00  )  00  TO  29 

5430-  Pin-PJUKI  /  NU) 

5440-  RP1  -  R  ( I )  1  Rdf  1  ) 

5450-  UMAX  -  TAU  *  SORT (  RP1  /  RM.1.  ) 

54 AO-  UMAX2  -  UMAX  »  UMAX 

5470-  If  ( R  ( I )  .  I..T  .  RSFEAR (  NU  )  )  GO  TO  20 

5480-  GO  TO  35 

5490-  20  BINGO- 1 . 

5500-  GUS-G ( UMAX ) 
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5500=  GUS=G  ( UMAX ) 

55!  '-  f  ? :  f ' .  nan  w,>: ) 

51.:-;  lii:  , a 

5550=  25  S  1 1  IGU  •  (  RSFEA  R  ( NU )  /R  Cl  )  /  .it  2 

:  55 AO  =  US  *  (  T AU/RM l  >  *  < SQR  I'  <  R  <  I  +  J  )  ##2-RS2  >  -SORT  (  R  ( I  )  **2~RS2  )  > 

5550-  HUS  =  1 1  (US) 

i  5500  =  GUS  ==  (.(115) 

5570=  28  F  (2.1  !  1  ,NU )  =  P1  Bi  ( RPl  1 1  AU  )  t*2/ ( A  .  *R  ( I  i  1 )  **2)  >fc<  (  ( ( F  ( 2 . 1 

:  5500  =  1  .ML!)  -  F  ( :L .  I  .MLJ)  )  /  SINSO)  t  ( G  (  TAIJ )  -  GUS  )  )  +  ( F  (  1 ,  I ,  MU )  -  PI  B 

5590==  2)  *  (GUS  -  GUJnAX)))  (■  RM1**2  /  (  2 .  *  (  R  ( If  1  )  *TAU )  «2 )  *  (  ( F  <  1 . 1 

5600==  3t  :L  .  NU  )  -  P I B )  *  ( 1 .  -  H  ( UMAX 2 )  )  -  (  ( F  (2,1. NU )  -F  <  1 » I » MU  )  )  / 

5610=  40INSN  )  *  <  H  (  TAU  )  -  HUS )  --  ( I"  (  I  » I ,  NU )  ~P  IB)  *(HUS  -  H(UMAX))) 

5620“  F  ( 2  » 1 1  1 .  NU )  =  ABS  <  F  (2.1 1 1 » NU  )  ) 

5630 =  GO  TO  30 

56 A 0“  29  F  ( 2  >  1  M. NU )  “ABC  <  8 .  *PI *  ( B  ( I » NU )  -B  ( 1  +  1 ,  NU )  )  /  ( 3 .  *  ( TAU IT A UP  )  )  ) 

5650==  30  CONTINUE 

5660=  X  ~  0. 

i  5670=  DO  35  NU  ==  l.MAXNU 

;  5680-=  35  X  =XfR  ( 2 )  «2*  ( F  (2.2.  MU )  -F  (1.2.  NU )  ) 

j  5690=  DEBT  ( .1. )  =3.  *  X  /  R<2>«3 

;  5700=  DO  50  I  =  2.100 

5710=  X  =  0 « 

5720=  DO  AO  NU  =  l.NAXNLJ 

5730=  40  X  =  X  +  <R<I)**2  *  (F  (2. 1 . NU) -F ( 1 . 1  .NU ) )  -R(I-fl)**2* 

!  57 A 0=  1  (F(2. 1  +  1  .HLJ)  -•  F(1.1  +  1.NU)  )  ) 

i  5750=  50  D  E  D  T  ( I )  =  ( -  X  %  3 .  /  ( R  ( I  i  1 )  *  *  3  -  R  ( I )  *  *  3 )  ) 

'  5760=  DO  60  1=1.100 

5770=  ER(I)  =  DEBT ( I ) 91  iT  i  ER(I) 

;  5780=  IF  ( ER  (I )  .  LT  .  6 . 122E--5 )  ER  ( I )  =6. 122E-5 

l  5790=  60  CONTINUE 

;  5800=  RETURN 

i  5810=  END 

i  5820=0 

;  5830=  FUNCTION  GUJ) 

‘  5840=  REAL  HMD El 

j  5850=  IF  (U. GT .225) GO  TO  10 

!  5860=  X=  EXP(-U)  *  ( 1  ,/U**2  -  1./  U) 

I  5870=  Y -  MMDEI (2.U. IER) 

j  5880  -  G=  X  ■{■  Y 

:  5890=  RETURN 

!  5900=  10  G "0 . 

|  5910=  RETURN 

i  5920=  END 

I  5930  C 

i  5940=  SUBROUT INF  PLANCK ( MAXNU . TEMP . B » WPTEMP . AVON!) ) 

5  5950=  D I  MENS  I  ON  TEMP  (  1 0 1 )  .  B  ( 1 0 1  .  J.  0.)  .  F  ( 9 ) »  AVGNU  (10) 

5960=  REAL  K 

5 9 7 0 C -  B - •  E R G /  <  S E C - •  C M **2>»  H = E R 6 - SE C  .  K = E R G / K E LVI N .  X  I'  - 1 1 H = A  V G N U  <  I  >  =  1 1 

5980  r.i  A  T  A  H .  K .  C  / 1 . 6  2 6 E  -  2  7 . 1 . 3  S 1 E  -16, 3  .  E 10/ 

5990=  A  ( U ,  L. )  =  II  *  U  /  (K  *  TEMP  ( l  )  ) 

6000=  PK(U.R)  =  2.  *  I  I  *  U**3  /  ( GY* 2  *  (EXP(R)  -•  1 )  ) 

6010=  X  =  UIP  LEMP  *  2.41800  /  FLOAT  ( MAXNU ) 

\  6020=  Y  =  WITF.MP  *  2. 41  BE  15 

6030-  llll  =  (X--Y)/n. 

6040  AVGNU <  I  )  =  ( X1Y ) /2  . 

6050=  IF ( MAXNU. EQ. J )  GO  TO  30 
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60 '.10  II  <  M A >: M II .  E 0 .1)  0 0  T 0  3 0 

6060  •  III)  *.I  I  -  2y|iAXNU 

60/0  -  5  AVGNU(I)  -  AOGNIK  1-1)  f  X 

6000  -  DO  20  NU  -  I.  >  MAXNIJ 

f  . [ ,  ■  r i  .  "/  y 

;  y  1  ;  >. 

6:1  I  0  -  l.i C)  20  I  ~  i  y  100 

6  J  20-  D  -  A  ( Z  y  I ) 

6 1  3  0  -  I  r  (1,1.  or.  2  2 5 .  )  F  <  1  )  -  0 

6140-  IF ( D . LE . 225 . )  F(l>  -  PK(ZyD) 

61110  -  D  -  A  <  Y  y  I ) 

6 1 60-  I F  <  H . GT .225.)  F ( 9 )  -  0. 

j  6 1 70-  I F  ( I! .  LEI .  225 .  )  F  <  9  )  -  PK  ( Y  r  D ) 

,  6180=  HO  10  J-  2  y  0 

I  6190=  Z  -  Z  +  Mil 

;  6200-  D  -  A ( Z  y I ) 

j  62 1 0-  I F  ( II .  GT  .225 .  )  F  ( J )  -  0  . 

i  6220-  I F  <  D . LE . 225 . )  F ( J )  -  PK ( Z  r D ) 

i  6220---  io  continue: 

•  6240-  20  B  ( I  y  NU )  -  HI  1/3.  *  ( F  <  1 )+  F  ( 9  >  +  4  .  *  ( F  ( 2  )  +F  <  4  )  +F  ( 6 )  +F  <  8 )  )  +2 .  * 

!  6250-  1  ( F  <  3 )  +F  ( 5 )  -IF  (?)))/( X  -Y ) 

i  6260-  GO  TO  50 

:  6270-  30  BO  40  I  -  ly.1.00 

:  6280-  B  ( I  y  1 )  -5 , 6686E-5*TEMP  <  I )  «4 

;  6290-  40  continue: 

;  63oo=  50  continue: 

•’  6310-  RETURN 

;  6320=  END 

^  6330= C 

•  6340-  FUNCTION  AMU  ( D  y  T  y  ANU  y  RI  IO  > 

6350-  REAL  MU (5) 

:  6360=  DIMENSION  TEMP (5) 

I  6370-  DATA  (TEMP (I ) y 1  =  1 »5) » ( MU ( J ) y  J-l y  5 ) /l . E3  y 1 . E4  y 1 . E5  y 1 . E6  y 1 . E7  y 1 .E- 1 

>  6380-  1  y  1 . 9 IE- 10  v  1 . 268E-8 1  J.  .  26E--6  r  1 . 27E-6/ 

I  6390-  IF  (T.LE.1.E3)  GO  TO  10 

;  6400-  IF  (T.GE.l.E?)  GO  TO  20 

•;  6410-  GO  TO  30 

;  6420-  10  AMU-7. 145E-"2*D/RI  iO 

i  6430-  RETURN 

j  6440-  20  AMU -6 . 25E-9*D/RI  10 

6450-  RETURN 

6460-  30  1=1 

6470-  IF ( T . GT ♦ 1 ♦ E4  >  1=2 

!  6480-  IF<T.GT, 1 ,E5)  1-3 

6  4  9 0  -  I F  <  T  .  G  T  .  1  .  E  6 )  I  -  4 

•  6500-  X-AI..OG ( I  .  /MU ( I )  )■!.(<  AI...OG (  TEMP ( I  !  1  )  )~ALOG(T)  )/<ALOG (TEMP <1  +  1  >  ) 

:  6510-  1  -ALGO  ( TEMP  <  I)  >  )  )  TAI...OG  ( MU  ( III  )  >  *  ( ALOG  (T ) -ALUG  UV.':  ;P  <  .1.  )  )  ) 

6 5  2 0  -  2 /  ( A  I  .  C)  0  <  T  E  Ml'  ■.  I  » .!.  )  1  A  L  0  G  (  I  E  M  P  ( I )  )  ) 

6530-  IF < X , GT . 225 )  X  225 . 

6540  -  AMIR  EXP  (  -X )  TB/EIiO 

6550-  RETURN 

6560  -  END 

6570  G 

6580  -  subruui  inf  gadget < maxnij y nr  y  wit  emp y  temp y  tstar j  y  n y li  y  ykt 

6590-  1  y  Ik  t  AMASS  y  R  ) 

6600-  I  I  1  M  L  N  G 1 0  N  A  r  i  A  8  G  (  J.  0 1  )  y  R  (101) 
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6600= 

n i mens ion  amass <  i o i )  * r <  i o  i.  > 

(Vi. 1.0 

ninnioi  no  rrupdoi  '•  .mioi  ).-mioi  >  •,  fr<  i  oi ' 

66 

c.lil-  l  !  I,-  ::  TINTS- UPTFMP  REV,. i  G/CM3  - - 

6:,  .,0: 

U 0  10  1  -  2,100 

6640= 

f  K(  1  >=6,  :l  7 IT. “5 

66l.»0“ 

CI(I>  =  3.8001: 6 

6660:- 

A  M  A  0  0  (  1.  )  =  0(1)  :i:  4  .  /  3 .  *  A  C  0  S  <  1  .  )  XT  R  < I  + 1  )  *  1 3 -  R  ( I )  *  *  3 ) 

6670“ 

10 

TEMP ( I )  =  300. 

6680= 

TSI  ART  =  1.E--7 

6690- 

READ*  ,  MAX  NO 

-  6700= 

PRINT*,  "  NUMBER  OF  FREQUENCY  GROUPS  ISJ  ’ ,MAXNU 

1  6710= 

C  A  L  L,  T  E  M  P  0  (  Y  K  T  ,  R  (  2 )  ,1  ,  El  ZERO,  E  R  2  E  R  0  ) 

6720= 

W  P  T  E  M  P = B  .  6 1 7  E  -  8  *  T 

6730= 

PRINT* , “ INITIAL  CELL  TEMP  =  " , T , *  KELVIN ,  OR  *  » UP TEMP, *  KEV " 

67-40= 

El ( 1 >=EIZERO 

:  6750  = 

TEMP ( 1 ) =T 

j  6760= 

ER < 1 > =ER7ER0 

6770= 

D(1 )= 0.33 86 

6780  = 

AMASS ( 1 ) =1 . E6 

6790= 

RETURN 

6800  = 

END 

6810=C 

■  6820= 

SUBROUTINE  TEMPO  <  YKT »R , T, D ,C) 

6830  = 

R  =  89. 

6840= 

A=YKT*4 . J 8E19 

6850  = 

TL0=1 . E3 

6860= 

THI---1  .E8 

6870  = 

10 

T=  ( THI  +  TL.0 )  /2 . 

.  6880= 

B=2.E14*T 

6890= 

C  =  5 , 6  6 8 6 E  ~ 7  *  T  *  *  4 

'  6900= 

b  =  c  •}■  n 

;  6910= 

IF(A.GT.B)  TLO=T 

;  6920= 

IF(A.LT.B)  TIT  I  =  T 

:  6930= 

I F ( ABS ( A-B ) . GE . 0 . 0001  * A )  GO  TO  10 

|  6940= 

C =7.5  5 8 E  -•  1 5  *  T  *  *  4 

:  6950= 

B=1 . 103E4*T*16.4 

!  6960= 

20 

RETURN 

j  6970= 

END 

•  6980=C 

{  6990= 

S (J If  R 0 IJTIN E  W R 1 7  E R  (  T E M P ,  P ,  D ,  E I ,  X ,  D E D T  ,  E R ,  F ) 

j  7000= 

DIMENSION  TEMP (101) , DEBT ( 101 ) ,EI ( 101) ,1(101) ,0(101) ,X(101) 

;  70io= 

I.i  I  MEN  S 1 0  N  E  R  ( 1 0 1 )  ,  F  ( 2 , 1 0 1  , 1 0 ) 

|  7020= 

PRINT  100 

i  7030= 

DO  10  I  =  1,20 

i  7040= 

J-l 

i  7050= 

Z=F (1,1, J ) 

:  7060= 

Y=F ( 2 , I , J ) 

■;  7070= 

10 

PR I NT  200 ,  I ,  X  ( I )  ,  TEMP  ( I )  ,  P  ( I )  ,  D  ( I )  ,  E I  ( 1 )  ,  DEBT  (  1 )  ,  ER  ( I )  ,  7. ,  Y 

!  7080= 

RETURN 

7090  = 

100 

FORMAT  < / / , 2  X , " I " , 10X, "X" , 10X, ’TEMP’ , 10X, ’PRESSURE* ,5X, *  DEN 

i  7100= 

ii ox*  "Ei" , i o x ,  “ heut ■  ,.i.ox,  *er"  ,.i.ox,  ’f<->  • , .1  o\ ,  "i  ;  • ) "  ,//> 

7110  = 

200 

FORMAT  (  2 X ,  1. 2 » 9  ( 2X » E .1. 2 . 5 )  ) 

7120= 

END 

71 30=C 

.  7140= 

FUNCTION  TEMPER (E I, T,V) 

.  7150= 

Y=28TAR<  T )  »  1. 
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7 j  y  ;::n (  i  >  t  i . 

7  i  •  ■  v-  :i  . 

7170-  if-  ( t  . c;  i ,  2 . 1:3 )  «=7. 

7 1 1;0“  TEMPER  I  .  1 69E-7XH 1  *  V/ ( 0* Y ) 

7190  RETURN 

7200-  END 

721 0--C 

7220 »  FUNCTIO M  Z S T  A R  ( T  ) 

7230-  IF(T.LE.2.E4>  ZSTAR=0  * 

7240*  I F ( T . GE *  1 . E6 )  Z STAR  =  7 . 

7250-  IF  Cl  .GT.2.E4. AND.T.LT. 1 .E6>  GO  TO  10 

7? AO*  RETURN 

7270-  10  ZSTAR-:  ( AL0G10 (  H-4.301  )  *4  . 

7280=  RETURN 

7290*  END 

7300  -C 

7310*  FUNCTION  H<U> 

7320=  IF ( U ♦ GT . 225 . )  GO  TO  10 

7330=  H---EXP  (  -1J )  *  <  1  .  -HJ ) 

7340=  RETURN 

7350=  10  H=0. 

7360  =  RETURN 

:  7370=  END 

:  7380=0 

-  7390=  FUNCTION  El A < P y D , RHO ) 

•  7400=0 -  GUESS  GAMMA  MINUS  ONE  - 

7410=  OKI =.20 

•  7420=  DO  30  1=1 y 10 

.  7430=0 -  COMPUTE  AN  INTERNAL  ENERGY  - - 

;  7440=  EI=P/<bM:l.*D> 

.  7450=0 -  COMPUTE  A  NEW  GAMMA  MINUS  ONE  - 

7460=  GM1  =  VGA  MM  A  (Ely  D  *  Id  10  )  - 1 . 

:  7470  =  0 - - -  COMPUTE  A  PRESSURE  - 

;  7480=  P2=6MJ*D*EI 

.  749o=C -  COMPARE  PRESSURES  - - - 

7500=  I F  <  ADS  ( P--P2  )  .  LE .  .001  *P )  GO  TO  40 

|  7510=0 - -  GAMMA  IS  CORRECT  WHEN  THE  PRESSURES  AGREE - 

•  7520=  30  CONTINUE 

;  7530=  PRINT** "I*"* I 

:  7540=  40  El A-ET 

7550=  RETURN 

7560=  END 
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